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STATEMENT  OF  PROBLEM 


The  ultimate  goal  of  this  research  is  to  insure  proper  design  of  molecular  sieve  oxygen 
generation  systems  for  the  U.  S.  Army’s  in-flight,  medivac,  and  Field  hospital  use. 
Specifically  the  research  involves  further  development  of  an  OBOGS  model  to  include  the 
effects  of  contaminants  in  the  feed  air.  This  OBOGS  model  can  be  used  to  optimize  and 
design  OBOGS  systems  with  respect  to  system  parameters  such  as  cycle  time  and  bed  and 
valve  dimensions. 


t. 


BACKGROUND 


The  pressure  swing  adsorption  process  (PSA)  and  its  utility  for  oxygen  generation  in 
military  applications  were  discussed  in  the  last  annual  report  (dated  September  13,  1986). 
PSA  can  also  be  used  to  generate  nitrogen  (using  zeolite  4A  in  the  packed  beds),  which  is 
useful  for  fuel  tank  inerting.  See  the  last  annual  report  for  a  discussion  of  the  PSA 
process. 

A  major  development  since  the  last  report  has  been  the  transfer  of  the  experimental 
work  from  USAARL  to  UT.  Thus,  much  effort  has  been  expended  to  develop 
experimental  techniques  that  will  allow  the  determination  of  unknown  model  parameters. 


I'.- 


' -■> .•> .•  >  >.> 


KK5 


TO 


39 


APPROACH  TO  THE  PROBLEM 

Chromatographic  experiments  will  be  done  to  estimate  model  parameters.  These 
experiments  should  yield  information  on  both  equilibrium  and  diffusion  phenomena.  A 
mass  spectrometer  has  been  ordered  and  should  be  delivered  at  the  end  of  this  month.  This 
device  will  potentially  allow  for  determination  of  ternary  equilibria,  where  the  use  of 
isotopically  labeled  components  may  be  necessary.  Once  the  experimental  arrangement  has 
been  tested  against  known  data,  work  can  begin  on  contaminants  with  unknown 
adsorptive  properties.  Experimental  apparatus  have  been  designed  for  both  single  column 
and  dual  column  experiments. 

A  new  equilibrium  model  will  be  implemented  to  allow  for  the  use  of  more  general 
isotherms  in  the  overall  model.  This  will  eliminate  the  assumptions  implicit  in  the  coupling 
scheme  presently  used.  The  new  model  can  predict  multicomponent  isotherms  from  pure 
component  data.  A  variation  of  this  model  contains  terms  for  "tuning"  the  isotherms 
generated  if  adequate  experimental  data  exist. 

In  addition,  the  mass  balance  equations  will  probably  be  modified  to  include  axial 
mass  dispersion.  Dispersion  may  be  important  at  the  zero  velocity  point  in  the  bed. 

A  new  computational  scheme  will  be  required  to  solve  the  new  system  of  equations. 
Adaptive  grid  finite  element  methods  are  being  considered. 
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RESULTS  AND  CONCLUSIONS 

Progress  has  been  made  in  two  areas:  experimental  and  computer  modeling.  The 
experimental  work  has  just  begun,  while  the  computer  work  is  a  continuation  and 
expansion  of  previous  efforts. 

The  present  computer  model  performs  well  for  the  oxygen-nitrogen  system.  This 
good  performance  is  the  result,  to  some  extent,  of  the  fact  that  these  gases  are  very  similar 
in  their  physical  properties.  The  system  is  equilibrium  dominated,  so  diffusional  effects  are 
minimized.  Similar  molecular  weights  and  diameters  further  "simplify"  the  calculational 
load.  Recent  modifications  fixed  some  compiler-specific  "bugs"  (only  discovered  when  the 
model  was  compiled  from  source  code  at  a  remote  location)  and  included  a  power  law 
correlation  for  the  isotherm  parameters  with  temperature  (see  Appendix  A.)  (1).  This 
temperature  correlation  extends  the  range  of  applicability  of  the  model.  However,  these 
recent  modifications  do  not  address  the  problems  that  can  be  expected  as  the  model  is 
applied  to  more  complex  systems.  The  model  is  being  rewritten  to  include  multi- 
component  feed  mixtures,  new  equilibrium  expressions  and  multiple  mass  transfer 
coefficients. 

The  number  of  components  in  the  feed  stream  is  presently  fixed  at  two  (specifically, 
oxygen-nitrogen).  This  is  to  be  expanded  to  three.  The  properties  of  the  third  component 
(contaminant)  are  likely  to  differ  greatly  from  those  of  oxygen  and  nitrogen.  Therefore, 
simplifications  possible  in  the  oxygen-nitrogen  system  (such  as  simple  isotherm  forms  and 
mass  transfer  expression)  will  not  apply  to  the  ternary  system. 

In  the  present  model,  oxygen  is  assumed  to  have  a  linear  isotherm  while  nitrogen  has 
a  simple  Langmuir-type  isotherm  (2).  The  isotherm  parameters  are  calculated  from  an 
analysis  that  is  specific  to  these  isotherm  forms.  The  isotherms  are  coupled  using  the 
concept  of  the  Ideal  Adsorbed  Solution  (3).  The  Ideal  Adsorbed  Solution  Theory  (IAST)  is 
valid  only  for  similar  molecules,  and  is  computationally  intractible  for  ternary  solutions. 

As  a  preliminary  step  to  a  general  ternary  model,  a  limited  ternary  model  has  been 
developed  (see  Appendix  B.).  In  its  present  form,  this  model  uses  uncoupled  linear 
isotherms  with  a  common  mass  transfer  coefficient.  The  isotherm  slopes  are  not  correlated 
with  temperature  This  model  represents  a  limiting  case  as  the  simplest  possible  ternary 
model.  Unfortunately,  it  has  not  been  tested  due  to  a  lack  of  experimental  data. 
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The  general  expansion  to  three  components  presents  problems  in  two  equilibrium 
areas:  the  form  of  the  pure  component  isotherm  for  the  third  component  and  the  method  of 
coupling  the  isotherms.  Because  the  contaminant  will  most  likely  be  present  at  low 
concentrations,  a  linear  isotherm  might  be  used.  However,  the  contaminant  will  be 
probably  be  quite  different  from  oxygen  and/or  nitrogen  so  IAST  w ill  not  apply.  Also, 
some  highly  polar  contaminants  may  exhibit  non-linear  behavior  even  at  low  concentrations 
due  to  localized  interactions  in  the  zeolite. 

A  literature  search  has  provided  a  model  that  appears  to  be  able  to  solve  the  pure 
component  isotherm  and  coupling  problems.  The  Vacancy  Solution  Model  (VSM) 

(4,5,6, 7)  includes  terms  for  both  adsorbate-adsorbate  and  adsorbate-adsorbent  non¬ 
ideality.  The  VSM  can  be  correlated  with  temperature  and  can  include  information  from 
available  binary  isotherms  to  improve  the  ternary  prediction.  Calculation  of  activity 
coefficients  may  be  time  consuming,  however,  in  the  VSM.  A  sensitivity  study  may  be 
done  to  determine  if  these  need  to  be  calculated  at  every  time/distance  step. 

A  preliminary  code  for  the  basic  VSM  (4,5)  has  been  written,  but  not  yet  fully  tested. 

In  the  present  model,  the  dynamics  of  mass  transfer  are  represented  in  the  classical 
form:  constant  mass  transfer  coefficient  with  linear  driving  force.  The  coefficient  is  the 
same  for  both  oxygen  and  nitrogen;  its  value  is  determined  empirically.  For  the  oxygen- 
nitrogen  system,  performance  is  relatively  insensitive  to  the  value  of  the  mass  transfer 
coefficient  because  the  system  is  equilibrium  dominated  (diffusion  is  rapid  for  both 
components).  The  system  equations  are  simplified  because  of  the  common  mass  transfer 
coefficient. 

Larger  contaminant  molecules  will  not  diffuse  as  rapidly  as  oxygen  and/or  nitrogen. 
Therefore,  at  the  very  least,  the  third  component  will  need  to  be  considered  separately.  The 
addition  of  this  equation  may  make  the  model  stiff  and  more  difficult  to  calculate.  This 
equation  will  require  a  value  for  the  mass  transfer  coefficient  for  the  third  component.  At 
present,  a  literature  search  is  underway  to  determine  the  best  approach  to  predicting  this 
coefficient,  as  well  as  the  form  of  the  equation.  It  is  most  likely  that  the  present  form  will 
be  retained  with  a  separate  expression  for  each  component.  However,  models  with  the 
mass  transfer  resistance  lumped  at  the  zeolite  crystal  surface  and  bidisperse  models  (with 
resistances  at  both  the  pellet  and  crystal  surfaces)  are  also  being  studied. 

The  inclusion  of  axial  dispersion,  along  with  the  changes  mentioned  above,  will  make 
a  new  numerical  scheme  necessary.  Finite  element  methods  are  being  considered.  Because 
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of  the  near  "shock"  phenomena  expected  in  the  concentration  profiles  through  the  beds,  a 
very  fine  grid  may  be  required  for  adequate  resolution.  Unfortunately,  fine  grids  are  very 
time  intensive  to  calculate.  Two  adaptive  grid  methods  are  under  review.  One  uses  preset 
"levels"  of  adaptation  to  cause  fine  grids  to  be  used  where  necessary;  more  coaise  grids  are 
used  elsewhere.  This  method  appears  to  suffer  from  the  fact  that  nodes  are  being  created 
and  lost,  which  may  complicate  its  implementation.  The  other  method  uses  a  constant 
number  of  nodes;  the  nodes  are  "moved"  to  form  a  tight  grid  where  it  is  required.  This 
method  suffers  from  problems  associated  with  the  equations  of  motion  for  the  nodes. 

The  main  structure  of  the  program  will  be  changed  to  provide  more  modularity.  This 
modularity  is  important  in  order  to  allow  for  flexibility  in  designing  an  OBOGS  with  the 
simulator.  As  three  bed  systems  and  systems  with  stratified  beds  (with  layers  of  different 
zeolite  types)  become  common,  the  model  will  be  required  to  simulate  them.  In  its  present 
form,  this  will  necessitate  a  major  rewrite  of  the  code  for  each  change.  A  more  modular 
approach  will  make  model  growth  more  convenient  and  efficient.  Specifically,  if 
temperature  effects  are  to  be  addressed  (i.e.,  a  non-isothermal  bed)  in  a  future 
implementation  of  the  program,  it  will  be  a  great  benefit  if  the  bed  equation  section  can  be 
easily  removed  and  replaced  instead  of  rewriting  the  entire  code. 


The  first  structural  change  will  be  the  movement  of  bed  structure  functions  from  the 
simulator  program  to  a  utility  program.  The  designed  bed  systems  w  ill  be  stored  as  data 
files  for  use  by  the  simulator.  This  will  remove  much  "overhead"  from  the  simulator  as 
well  as  reduce  the  number  of  times  the  program  needs  to  be  recompiled  when  making  bed 
parameter  changes.  A  graphics  library  (for  the  DEC  Regis  VT  terminals)  as  well  as  a 
possible  structure  for  utility  data  files  have  been  completed. 
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Isotherm  information  for  many  pure  substances  on  zeolites  is  not  readily  available. 
Diffusion  information  is  even  more  rare.  Air  separation  data  on  zeolites  5A,  4A,  and  1 3X 
are  available  in  the  literature  because  the  process  is  important  and  common.  However, 
information  on  other  substances  is  either  completely  lacking,  at  inappropriate  conditions 
(for  instance,  in  catalytic  cracking  processes  which  operate  at  elevated  temperatures  and 
pressures),  or  "proprietary.” 


Therefore,  experimental  work  must  be  done  to  determine  values  for  necessary  model 
parameters.  Chromatographic  techniques  hold  out  the  promise  of  obtaining  both  isotherm 
and  diffusion  information  relatively  quickly  and  simply  (8,9).  Isotherm  parameters  can  be 
obtained  by  measuring  pulse  retention  volumes  at  low  flowrates.  Diffusion  data  can  be 
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obtained  from  a  moments  analysis  on  the  outlet  peak  at  higher  flowrates.  The  moments 
analysis  usually  assumes  a  symmetric  peak  (for  a  pulse  experiment)  and  linear  isotherms, 
however;  its  applicability  to  the  OBOGS  problem  will  be  unknown  until  some  preliminary 
experiments  have  been  done. 

Most  of  the  chromatographic  approaches  to  determination  of  isotherm  parameters 
usually  depend  on  an  assumption  of  the  linear  isotherm  for  the  components  at  a  given 
composition.  For  small  enough  concentration  pulses,  this  may  be  a  reasonable 
assumption.  The  assumption  of  a  linear  isotherm  for  a  concentration  pulse  is  best  when  an 
isotope  of  one  of  the  components  is  used  as  the  pulse  (9).  The  pulse  can  then  be  large 
(because  the  bulk  concenetration  will  not  change,  no  matter  what  the  pulse  size)  and, 
consequently,  more  easily  detected.  The  slope  of  the  isotherms  at  each  composition  (which 
are  the  data  obtained  from  the  chromatographic  analysis)  must  be  integrated  to  obtain  the 
isotherm  itself.  Then,  the  isotherm  parameters  can  be  estimated  from  the  isotherm  (for  the 
VSM  model  (without  temperature  correlation)  there  will  be  three  parameters  to  estimate). 
The  fitting  procedure  used  will  have  a  critical  effect  on  the  quality  of  the  parameters 
estimated.  Traditional  least  squares  analysis  generally  assumes  all  of  the  error  is  in  the 
dependent  variable.  In  our  work,  however,  the  independent  variable  (partial  pressure)  is 
also  measured  and  subject  to  error.  A  recent  paper  (10)  has  noted  that  while  the  least 
squares  method  will  fit  the  pure  component  isotherms  well,  the  fact  that  the  error  in  the 
independent  variable  was  ignored  will  cause  significant  deviations  when  one  attempts  to 
predict  binary  (or  higher)  isotherms  from  these  pure  component  parameters.  Therefore,  the 
Frror-in- Variables  Method  (EVM)  is  being  considered  for  use  in  parameter  estimation  (11). 
This  method  accounts  for  error  in  both  measurements  through  the  use  of  covariance 
matrices  based  on  many  experiments.  Construction  of  these  matrices  would  be  quite  time 
consuming.  However,  it  appears  as  though  published  covariance  matrices  based  on  similar 
experiments  may  be  applicable  to  our  work  (10). 

Some  preliminary  work  on  the  isotherm  analysis  has  been  done  already,  using 
available  equipment  at  Brooks  Air  Force  Base.  These  experiments  were  done  in  an  attempt 
to  duplicate  previous  work  (8)  in  order  to  validate  the  method.  Due  to  equipment 
limitations,  only  the  binary  ^  ./gen-nitrogen  isotherm  was  determined.  Unfortunately, 
during  the  analysis  of  these  experiments,  errors  were  found  in  the  mathematical  derivations 
in  (8)  which  led  to  non-physical  results.  The  equations  were  solved  again  and  physically 
realizable  isotherms  were  obtained  (the  corrected  equations  are  in  Appendix  C).  However, 
as  the  analysis  in  (8)  was  faulty,  it  was  decided  not  to  compare  our  results  with  those  in  the 
paper.  A  comparison  to  volumetric  data  on  this  mixture  will  be  done  (it  would  be  better  to 


compare  to  data  on  the  mixture  at  constant  total  pressure,  but  binary  data  of  this  type  is 
rare)  to  check  on  the  validity  of  the  method. 

Single  and  double  bed  apparatus  have  been  designed  for  the  anticipated  experiments. 
A  mass  spectrometer  has  been  ordered  and  is  expected  to  arrive  at  the  end  of  this  month. 
Experimental  work  can  then  begin  here  at  UT. 
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APPENDICES 


A.  Isotherm  Temperature  Correlation  Equations 

B.  Simple  Ternary  Feed  OBOGS  Model 

C.  Corrected  Equations  from  van  der  Vlist 


APPENDIX  A. 


The  temperature  correlation  of  the  isotherm  parameters  was  entirely  empirical  work 
done  by  Walshak  (1).  Subroutines  CAP1,  INITIAL,  and  MENU  need  to  be  modified  to 
implement  the  temperature  correlation;  they  are  included  in  this  appendix.  The  temperature 
correlations  are: 

For  oxygen: 


Ka=  .1423275  +  1.83744e-3  *  (T)  +  7.2443e-6  *  (T)2 
For  nitrogen: 


Kb=  .03654573  +  9.5775e-4  *  (T)  +  1.32le-5  *(T)‘ 
B  =  3.059125 


T  is  in  degrees  Celsius 


SUBROUTINE  CAPl 
REAL  KA,KB 

COMMON/CAP /KA , KB , B , D , TEMPI 
COMMON/VOID/E 

KA=. 1423276+ . 00183744-TEMP1+ . 0000072443-TEMPI  * -2 
KB= . 03664673+ . 0009677S-TEMP1+ . 00001321 -TEMPI - *2 
B=3 .0691260816 

WRITE (7, 660) TEMPI, B,D,E,KA, KB 
WRITE (6 , 660) TEMP 1 , B , D , E , KA , KB 
660  FORMAT (’  THE  CURRENT  BED  PARAMETERS  ARE:’,/, 

1’  12.’, T6,’  TEMPERATURE  =  ' , T30 , F8 . 3 , T40 , ’ DEGREES  C’,/, 

2’  13.’,T6,’  B  =  ’ ,T30,F8.4,/, 

3’  14 .  ' , T6 , ’  DIFFUSION  COEFFICIENT  =  * ,T30,F8 . 2 ,T40 , ’ 1/SEC  * , / , 

4’  16.*, T6,’  VOID  FRACTION  *  ' ,T30,F8 . 3 , / , 

6  T6 , ’  AT  THIS  TEMPERATURE,  KA  *  * ,T30,F8.4,T4O, 

8  ’ KGMOLES  02  ADSORBED/KGMOLE  GAS’,/, 

7  T6 , 21X, *  KB  =  ’ ,T30,F8. 4, T40, 'KGMOLES  N2  ADSORBED/KGMOLE  GAS') 
RETURN 

END 


SUBROUTINE  INITIAL  (Cl  ,C2  ,Nl ,  N2  P) 

REAL  L , MA , MB , NAOl , NA02 , KA , KB , N1 , N2 , NBOl , NB02 , LIN 
DIMENSION  N1 (300,2) ,N2(300,2) ,X(3) ,P(2) 

COMMON/BEDV/LUMPS , DZ , L , Dl , D2 , MA , MB , RT , AI , AO , AMW01 , AMWL 1 , AMW02 
1  ,AMWL2,WBR 

COMMON/CAP /KA , KB , B , D , TEMP 1 
COMMON/VOID/E 

COMMON/GEOM/DBYIN , DVSIN , DVOIN , LIN , DIAOIN , DIAIIN 

. . . . . . . . . . 

c 

C  RT  =  GAS  CONSTANT • TEMP  (N-M/KGMOLE) 

C  AITOT, AOTOT  =  INLET  AND  OUTLET  CONCENTRIC  BED  AREAS  (M**2) 

C  MA.MB  =  KGMOLECULAR  WEIGHTS  (KG/KCMOLE) 

C  O.L  =  SUBSCRIPTS  FOR  Z*0  AND  Z=L 

C  AMW  =  AVERAGE  MOLECULAR  WEIGHT  (KG/KGMOLE) 

C  1,2  =  SUBSCRIPTS  FOR  BED  1  AND  BED  2 

C  E  =  VOID  FRACTION 

C  TEMPI  =  TEMPERATURE  (DEGREES  C) 

. . . . 

RT=8314. * (TEMPI +273. 16) 

AIT0T-3. 14160* (DIAIIN* .0264) **2/4. 

A0T0T=3 . 14160* (DIAOIN* *2-DLAIIN*  *2) * . 0264**2/4 . 

L= . 0264*LIN 

PRINT  BED  LENGTH  (M)  =  1  ,L 
PRINT  INLET  BED  AREA  (M**2)  =  AITOT 
PRINT  *,'  OUTLET  BED  AREA  (M*«2)  =  ’ , AOTOT 
MA=32. 

PRINT  VOID  FRACTION  =  ’ ,E 

MB=28 . 

AI=AIT0T*E 

A0=A0T0T*E 

AMW01=MA»X(1) +MB* (1 . -X(l) ) 

AMWL1=MA*X (2) +MB* (1 . -X (2) ) 

AMW02=AMWL 1 

AMWL2=MA*X (3) +MB* (1 . -X(3) ) 

Cl=P ( 1) *6 . 806E3/RT 
C2=P (2) *8 . 806E3/RT 
CA1=C1*X(1) 

CA2=C2*X (3) 

CB1=C1-CA1 

CB2=C2-CA2 

NBOi= (CBl/KB) / (1 . +CB1/ (KB-B) ) 


NAO 1 =CA1 /KA 

NB02=  (CB2/KB)  /  (1 .  +CB2/  (KB*B)  ) 

NA02=CA2/KA 

DO  1  1=1, LUMPS 

Nl  (1, 1)=NA01+NB01 

Nl (I , 2) =NA01 

N2  (1 , 1) =NA02+NB02 

N2 (I , 2) =NA02 

DZ=L /FLOAT (LUMPS- 1 ) 

RETURN 

END 
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SUBROUTINE  MENU 

REAL  LIN ,  KA  , KB , LINl , KA1 , KBl 

DIMENSION  TI (11) , NY (11) 

COMMON/SYST/DT , PSUP , POUT , TF , WBRL , VS , VO , IIN , TCYC 
COMMON/GEOM/DBYIN , DVSIN , DVOIN , LIN , DIAOIN , DIAIIN 
COMMON/BEDV/LUMPS , DZ , L , D1 , D2 ,MA ,MB , RT , AI , AO , AMWOl , AMWLl ,AMW02, 
AMWL2 , WBR 
COMMON/VOID/E 
COMMON/CAP /KA , KB , B , D .TEMPI 
COMMON/NFINISH/NEXIT 

COMMON/ STEP /NSTEP , TSTEP1 , TSTEP2 , BRSTEP 
COMMON/NPRINT/NOUTl , TI , NY 
NEXIT=0 
NSTEP=0 
lO  CALL  SYSTl 
CALL  GEOM1 
CALL  CAPl 
WRITE (7, 107) 

WRITE (8, 107) 

107  FORMAT (’  ENTER  CORRESPONDING  PARAMETER  #  TO  CHANGE  PAR  VMETER ’ , / , 
1*  ENTER  "O"  TO  RUN  SIMULATION,’, 

2'  *68"  TO  STEP  THE  BREATHING  FLOWRATE’,/, 

4’  ENTER  "88"  TO  EXIT  PROGRAM’) 

N=0 

READ  *  ,  N 

IF (N. Eft.  O)  GOTO  800 
IF(N.EQ.SB)  GOTO  600 
IF (N . EQ . 88)  GOTO  800 
EPSI= . OOOOOOl 
WRITE (7, 120) 

WRITE (8, 120) 

120  FORMAT (/// ,  ’ . . . - . . 

1  //,’  ENTER  THE  PARAMETER  AND  PRESS  "RETURN"’,/, 

2/ /T40, ’CURRENT  VALUE’,/) 

IF (N . LT . 8) GOTO  lOO 
IF (N . LT . 12) GOTO  200 
IF (N . LE . 1 8) GOTO  300 
GOTO  800 
lOO  CONTINUE 

IF (N-2) 168 ,166, 166 
166  IF (N-3) 184 ,171,171 

171  IF (N-4) 170, 178, 174 
168  WRITE (8, 180)  PSUP 

WRITE (7, 180)  PSUP 

180  FORMAT ( ’  SUPPLY  PRESSURE  (PSIA) ’ ,T40 , F8 . 4 , /) 

READ  - .PSUPl 

IF (PSUP1 .LT . ESPI)  GOTO  lO 

PSUP=PSUP1 

GOTO  lO 

184  WRITE (8, 188)  POUT 

188  FORMAT ( ’  OUTLET  PRESSURE  (PSIA) ’ ,T40 ,F8 . 4 , /) 

READ  * , POUTl 

IF (POUTl . LT . EPSI)  GOTO  lO 

P0UT=P0UT1 

GOTO  lO 

170  WRITE (8, 172)  TF 

172  FORMAT ( ’  FINAL  OBSERVATION  TIME  (SEC) ’ ,T40 , F8 . 4 , /) 

READ  * ,  TF1 

IF(TFl .LT.EPSI)  GOTO  lO 

TF=TFl 

GOTO  lO 

176  WRITE (8, 177)  TCYC 

177  FORMAT (’  CYCLE  TIME  (SEC) ’ , T40 , F8 . 4 , /) 


READ  * , TCYCl 

IP (TCYC1 . LT . EPSI)  GOTO  lO 

TCYC=TCYC1 

GOTO  lO 

174  WRITE (6, 176)  WBRL 

176  FORMAT (’  BREATHING  MASS  FLOWRATE  (STD  LIT/MIN) * , T40 , F8 . 4 , /) 
READ  - , WBRLl 

IF (WBRLl . LT . EPSI)  GOTO  178 
WBRL = WBRLl 
GOTO  lO 
178  GO  TO  lO 
200  CONTINUE 

IF (N-7)  200,216,206 

206  IF(N-O)  226,238,207 

207  IF(N-IO)  236,246,266 
200  WRITE (6, 210)  DBYIN 

210  FORMAT ( *  BY-PASS  VALVE  DIAMETER  (IN) ’ ,T40,F8.4,/) 

READ  -,DBYIN1 

IF (DBYIN1.LT.EPS!)  GOTO  lO 

DBYIN=DBYIN1 

GO  TO  lO 

216  WRITE (6, 220)  DVSIN 

220  FORMAT ('  SUPPLY  VALVE  DIAMETER  (IN) * ,T40, F8 . 4 , /) 

READ  » ,DVSINl 

IF (DVSIN1.LT. EPSI)  GOTO  lO 

DVSIN=DVSIN1 

GOTO  lO 

226  WRITE (6, 230)  DVOIN 

230  FORMAT ('  OUTLET  VALVE  DIAMETER  (IN) * ,T40,F8.4, /) 

READ  • , DVOINl 

IF(DVOINl.LT.EPSI)  GOTO  10 

DVOIN=DVOINl 

GOTO  lO 

336  WRITE (6, 240)  LIN 

240  FORMAT (’  BED  LENGTH  (IN) ’ ,T40 ,F8 . 4 , /) 

READ  • , LINl 

IF (LINl . LT . EPSI)  GOTO  lO 

LIN=LINl 

GOTO  lO 

246  WRITE (6, 260)  DIAOIN 

260  FORMAT (’  OUTER  BED  DIAMETER  (IN) ’ ,T40 ,F8 . 4 , /) 

READ  • , D I AO INI 
IF (DIA0IN1.LT. EPSI)  GOTO  lO 
DIAOIN=DIAOINl 
GOTO  lO 

266  WRITE (8, 260)  DIAIIN 

260  FORMAT ( ’  INNER  BED  DIAMETER  (IN) ' ,T40 , F8 . 4 , /) 

READ  * , DIAIIN1 
IF (DIAIIN1.LT. EPSI)  GOTO  lO 
DIAIIN=DIAIIN1 
GOTO  lO 
300  CONTINUE 

IF (N- 13) 300 , 326 , 306 
306  IF(N-16)336, 346,10 
300  WRITE (6, 310)  TEMPI 

310  FORMAT ( ’  TEMPERATURE  (C) ’ , T40 , F8 . 4 , /) 

READ  - , TEMPO 

IF (TEMP2.LT. EPSI)  GOTO  10 
TEMP 1= TEMP 2 
GOTO  lO 

326  WRITE (8 , 330)  B 

330  FORMAT (’  B ’ , T40 , F8 . 4 , /) 

READ  *  ,  Bl 


IF(Bl.LT.EPSI)  GOTO  lO 

B=Bl 

GOTO  XO 

336  WRITE (0, 340)  D 

340  FORMAT (’  DIFFUSION  COEFFICIENT  (>100 . ) ’ ,T40, F8 . 4 , /) 

READ  * ,  Dl 

IF(Dl .LT.EPSI)  GOTO  lO 

D=D1 

GOTO  lO 

346  WRITE (6, 360)  E 

360  FORMAT (*  VOID  FRACTION  «1 . ) ’ ,T40,F8 .4, /) 

READ  * ,El 

IF (El . LT . EPSI)  GOTO  lO 

E=El 

GOTO  lO 

600  WRITE (8 ,6X0)  TF 
WRITE (7, 6X0)  TF 

610  FORMAT (////  . . ’//, 

1 ’  ENTER  THE  TIME  AT  WHICH  THE  STEP  CHANGE ’ , / 

3’  IN  BREATHING  FLOWRATE  IS  TO  OCCUR’,//, 

3*  THE  FINAL  OBSERVATION  TIME  IS  CURRENTLY  *,F7.3,’  SECONDS’,///) 
NSTEP*1 
READ  • , TSTEP1 
616  WRITE (8 , 618)  TF 
WRITE (7, 618)  TF 

618  FORMAT (////  ’ . - . - . ’//, 

1*  ENTER  THE  TIME  AT  WHICH  THE  STEP  CHANGE’,/ 

3’  IN  BREATHING  FLOWRATE  IS  TO  END’,//, 

3’  THE  FINAL  OBSERVATION  TIME  IS  CURRENTLY  ’,F7.3,’  SECONDS’,///) 
READ  • , TSTEP3 

WRITE (7, 630)  TSTEPl.TSTEP3.WBRL 
WRITE (8, 630)  TSTEPl , TSTEP3 , WBRL 

630  FORMAT  (/  /  ’  . . * . ••••*• . ’//, 

1 ’  ENTER  THE  NEW  BREATHING  FLOWRATE  OCCURING  AS ’ , / , 

3’  A  STEP  CHANGE  FROM  T=  ’.F7.3,’  TO  T=  ’.F7.3,’  SECONDS’,//, 

3'  THE  CURRENT  BREATHING  FLOWRATE  IS’,F6.3,’  STD  LIT/MIN ’ , ///) 

READ  * , BRSTEP 

WRITE (7, 630)  BRSTEP , TSTEPl , TSTEP3 
WRITE (8, 630)  BRSTEP , TSTEP 1 , TSTEP3 
630  FORMAT (//// , 

1’  . . . ’,////, 

3’  THE  BREATHING  FLOWRATE  OF  ’,F8.3,’  STD  LIT/MIN’,/ 

3’  WILL  BE  INPUT  AS  A  STEP  AT  T«  ',F7,3,  ’  SECONDS’/, 

3’  AND  WILL  END  AT  T*  ’,F7.3,  ’  SECONDS'////, 

4’  . ’////) 

GOTO  lO 
800  CONTINUE 

WRITE (7, 810) 

WRITE (8 , 810) 

810  FORMATC//////,  ’ . ’,/////, 

. ’  THE  FOLLOWING  PLOTS  ARE  AVAILABLE  TO  YOU  ON  THE  TERMINAL’,///, 

.’  1  MOLE  FRACTION  OF  OXYGEN  VS . TIME  (DATA  FILE=MOLEFRAC . DAT) ’ , / , 
.’  3  INLET  MASS  FLOWRATE  VS.  TIME  (DATA  FILE=FLOWRATE . DAT)  ’ ,  /  , 

.  ’  3  OUTLET  MASS  FLOWRATE  VS .  TIME  (DATA  FILE=FLOWRATE . DAT) ’ , /  , 
.’4  A  DYNAMIC  SIMULATION  FOR  ONE  CYCLE  (OXYGEN  VS.  TIME,’,/, 
DATA  FILE=PROFILE . DAT) ’ , 

•////, 

. ’  PRESS  THE  CORRESPONDING  #  TO  HAVE  THE  OUTPUT  SHOWN  AS  THE’,/, 

.’  SIMULATION  TRANSPIRES’ ,/////) 

READ (6 , 630) NOUTl 
830  FORMAT (II) 

IF (NOUTl .EQ. 4)  GOTO  760 
GO  TO  400 


760  DTI=1.0 

DO  J=l,10 
TI(J)=TI(J-1)+DTI 
NY  (  J)  =IFIX  (TI  (  J)  /DT) 

END  DO 

TI(11)=TCYC 

NY (11) =IFIX(TI (11) /DT) 

GO  TO  400 
400  WRITE(e,410) 

WRITE (7, 410) 

410  FORMAT (//// , ’  YOU  HAVE  THE  FOLLOWING  OPTIONS  :  * , 

•/////, TIO, *1.  ENTER  "O"  TO  RUN  SIMULATION.',/, 

. TlO ,  ’  2 .  ENTER  "10"  TO  CHANGE  THE  PARAMETERS.’,/, 

•TlO, ’3.  ENTER  "00"  TO  EXIT  PROGRAM.’,/////, 

. ’  ENTER  THE  CORRESPOrjDlNG  #  AND  PRESS  "RETURN"’,/////) 
read  *,  Nee 

IF(Nee.EQ.O)  GOTO  820 
IF(N60.EQ. 10)  GOTO  lO 
IF(NOO.EQ.OO)  GOTO  OOO 
800  CONTINUE 

WRITE (7, 810) 

WRITE(8 , 810) 

810  FORMAT (/////,  ’ . ’,/////, 

1 ’  THE  PARAMETER  VALUE  WAS  NOT  INPUT  CORRECTLY ’ , / 

2'  TRY  AGAIN’,/////,  ’ . ',////) 

GOTO  lO 

820  WRITE (8, 826) 

WRITE <7, 826) 

826  FORMAT (////, •  DO  YOU  WISH  TO  DOUBLE  THE  NUMBER  OF  SPACE  LUMPS  IN’/, 
1 ’  THE  SIMULATION  FOR  BETTER  ACCURACY  ?  (  CURRENT  NUMBER  USED  ’ / , 

2’  IS  60  )  (Y/N) ’,/////) 

READ  (6,830)  KEYIN 
830  FORMAT  (  A  ) 

IF  (KEYIN . EQ . YES)  THEN 

LUMPS  =  lOl 

ELSE 

LUMPS  *  61 
END  IF 
GO  TO  lOOO 
OOO  CONTINUE 
NEXIT=1 
lOOO  RETURN 
END 
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APPENDIX  B. 


This  appendix  contains  the  simple  ternary  feed  OBOGS  model  computer  code.  The 
model  uses  simple  linear  uncoupled  isotherms  for  all  components  with  a  common  mass 
transfer  coefficient 


PROGRAM  MSEV 


OBOGS  WITH  CONTAMINANT 
A  =  COMPONENT  A  (OXYGEN) 

B  =  COMPONENT  B  (CONTAMINANT) 

C  =  COMPONENT  C  (NITROGEN) 

REAL  LIN , KA , KB , KC 
REAL  L50 
BYTE  CONTINUE 
BYTE  LOWC 
BYTE  KEYIN' 

BYTE  NULL 
BYTE  YES 
BYTE  NO 

CHARACTER *84  TEMP 

DIMENSION  W (3) ,XA(3) ,XB(3) ,CVA(3) ,PS(2) ,P(2) 

DIMENSION  WV (502 , 2) , TV (500) , XV (502 , 3) 

DIMENSION  DL (600) , CA1 (600) , CA2 (600) , CBl (600) , CB2 (600) , 

1  XA1 (502, 11) , XA2(502, 11) ,XBl(502,ll) ,XB2(502,11) 

DIMENSION  IXABSl (601,11) ,IYABSl (601,11) ,IXABS2(601 , 11) , 

1  IYABS2 (601 , 11) 

DIMENSION  1X1(601) ,IY1 (601) , 1X2 (501) , IY2 (501) 

DIMENSION  JX (6,11) , JY(8,11) , JXERSH (8 , 1 1) , JYERSH(6, 11) , 

1  JXERSL(8,11) , JYERSL(8,11) 

DIMENSION  TI(8),  NY(8) 

DIMENSION  XCA1 (502, 11) , XCA2 (602 , 1 1) , XCBl (502 , 1 1) , 

1  XCB2(502,11) , XBR02(600) 

dimension  rangt(4)  ,xaxia(10)  ,  yaxia(lO)  ,ruigi2(4)  ,xaxia2(10)  , 
1  yaxia2(10) ,xajcia3(10) ,yajcia3(10) ,ruig*3(4) , 

1  x*xia60(10)  ,  yajciaSO  (10)  ,rang«60(4) 

DIMENSION  LSO (402) ,T60(400) , ERR (400) 

COMMON/TIME/T 

COMMON/SYST/DT , PSUP , POUT , TK , WBRL , VS , VO , IIN , TCYC 

COMMON/ GEOM/DBYIN , DVSIN , DVOIN , LIN , DIAOIN , DIAIIN 

COMMON/CAP/KA , KB , KC , BB , BC , D 

COMMON/VOID/E 

COMMON/NPINISH/NEXIT 

COMMON/STEP/NSTEP , TSTBPl , TSTEP2 , BRSTEP 
COMMON/NPRINT/NOUT1 , TI , NY 
COMMON/BEDV/LUMPS , DZ 
COMMON/HILOBED/NBED , NB 
DATA  CONTINUE/87/ 

DATA  LOWC/OO/ 

DATA  NULL/O/ 

DATA  YES/80/ 

DATA  NO/78/ 

CALL  MENU 

IF  (  LUMPS  ,E<J.  lOl  )  THEN 

IINC  =  3 

ELSE 

IINC  =  8 

END  IF 

IF  (  NEXIT  .EQ.  1  )  GOTO  COO 
IF  (  NSTEP  .EQ.  1  )  GOTO  20 
TSTEPl  =  lOOOOOOO. 

TSTEP2  *  lOOOOOOO. 

BRSTEP  -  WBRL 
20  CONTINUE 

DLIN  =  LIN  /  (LUMPS-1) 

DO  I  =  1, LUMPS 
DL (I)  =  (1-1)  *  DLIN 
END  DO 


I 


K*i 


OPEN  (UNIT= 1 ,  FILE= ’ MOLEFRAC  * , STATUS=  *  NEW ’ ) 

OPEN  (UNIT=2 ,  FILE= ’ FLOWRATE  1 , STATUS=  *  NEW ’ ) 

OPEN  (UNIT=3,  STATUS= ' SCRATCH ’ ) 

OPEN  (UNIT=4 ,  FILE= ’ PROFILE ’ , STATUS=  *  NEW  * ) 

OPEN  (UNIT=8,  STATUS= ’ SCRATCH ’ ) 

OPEN  (UNIT=9,  FILE= ’ REGIS ’ ,STATUS= ’NEW’) 

OPEN  (UNIT= lO , FILE= ’ CONTA ’ , STATUS= ’ NEW ’ ) 

WRITE (6 , 7) 

7  FORMAT (// , 

1  ' . ’) 

PRINT  * 

PRINT  *  ,’  OBOGS  SIMULATION’ 

PRINT  * 

IIN  =  O 
ND  =  O 
NSO  =  O 

CDBY  =  .056  /  DBYIN 

IF  (  CDBY  .GT.  1.  )  CDBY  =  1. 

IF  (  CDBY  .LT.  .6  )  CDBY  =  .6 

DBYM  =  DBYIN  *  .0254 

PRINT  BYPASS  VALVE  DISCHARGE  COEFFICIENT  =  ’ , CDBY 

BYVA  =  3. 14168-DBYM**2  /  4. 

VS  =  .8  *  (  .02S4-DVSIN  )**2  *  3.14168  /  4. 

VO  =  .8  «  (  .O264-DV0IN  )**2  *  3-14158  /  4. 

PRINT  AREA  OF  BY-PASS  (M**2)  =  * ,BYVA 

PRINT  *,’  CD* AREA  SUPPLY  VALVE  (M**2)  =  ’  ,VS 
PRINT  CD -AREA  OUTLET  VALVE  (U--2)  =  * ,V0 

PRINT  *,’  SUPPLY  PRESSURE (PSI) = ’ , PSUP , 

’OUTLET  PRESSURE (PSI) POUT 
CVA (2)  =  CDBY* BYVA 
T  =  0.0 
DT  =  .Ol 
DATAP  =  200. 

INC  =  IF  IX  (  TF  /  (DATAP  *DT)  ) 

NT  =  O 
IM  =  O 
P(l)  =  14.6 
P (2)  =  14.6 
XA(2)  =  .20 
XB (2)  «.  O. 

PRINT  «,’  BREATHING  FLOW  (STD  LIT/MIN)  =  * , WBRL 
XABR  =  . 20 
XBBR  =  O. 

IF (NOUTl -1 ) 36,35, 38 
36  CLOSE (UNIT=2 , STATUS= ’ DELETE ’ ) 

CLOSE (UNIT=4 , STATUS= ’ DELETE ’ ) 

CLOSE (UNIT=8 , STATUS= ’ DELETE ’ ) 

CLOSE (UNIT=10,STATUS= ’DELETE’) 

WRITE(1 ,67) 

67  FORMAT (’  OXYGEN  MOLE  FRACTION 24X ,  ’TIME  ’ ) 

GO  TO  1 

30  IF (NOUTl . EQ . 4)  GOTO  37 
IF (NOUTl . EQ . 5)  GOTO  170 
CLOSE (UNIT=1 ,STATUS=’ DELETE’) 

CLOSE (UNIT=4 , STATUS= ’ DELETE ’ ) 

CLOSE (UNIT=9 , STATUS= ’ DELETE ’ ) 

CLOSE (UNIT= 1 O , STATUS= ’ DELETE ’ ) 

WRITE (2 , 68) 

68  FORMAT ( ’  INLET  MASS  FLOWRATE ’ , T26 , •  OUTLET  MASS  FLOWRATE ’ , 
1’  T66 , TIME ’ ) 

GO  TO  1 

170  CLOSE (UNIT= 1 ,STATUS=’ DELETE’) 

CLOSE (UNIT=2 , STATUS= ’ DELETE ’ ) 


.•> 


37 


1 


77 

160 


CLOSE (UNIT=4 , STATUS= ’ DELETE  * ) 

CALL  FRAME 1 
CO  TO  1 

CLOSE (UNIT= 1 , STATUS= ' DELETE ’ ) 

CLOSE (UNIT=2 , STATUS= ’ DELETE  * ) 

CLOSE (UNIT=10 , STATUS = ’ DELETE ’ ) 

CALL  FRAME 1 
CONTINUE 
XA ( 1 )  =  .21 
XA  (3)  *  .21 
XB (1)  *  .21 
XB  (3)  =  .21 

CALL  TIMEF  (  CVA,PH,PS,T  ) 

CALL  BEDS  (  P,  W,  XA,  XB,  XABR ,  XBBR ,  CVA,  PS, 
CA1,  CA2,  CBl ,  CB2,  Cl,  C2  ) 

IF  (  NOUTl  .EQ.  4  .OR.  NOUTl  .  EQ.  6  )  GO  TO  77 
GO  TO  78 
DO  lOOO  K=l,8 

IF  (  NT  . EQ.  NY (K)  )  GO  TO  150 

GO  TO  lOOO 

CONTINUE 

IF  (  K  .GT.  1  )  THEN 

CALL  INITIALl  (  IYl(l),  IY2(1),  NB  ) 

DO  J=2, LUMPS 

DLIN  =  LIN  /  (LUMPS- 1) 

DL ( J)  =  (J-l)  •  DLIN 

IF  (  NB  .EQ.  O  )  THEN 

CALL  DRAWF  (  1X1 ( J) ,  IYl ( J) ,  1  ) 

ELSE 

CALL  DRAWF  (  1X2 (J),  IY3(J) ,  1  ) 

END  IF 
END  DO 

CALL  INITIAL2  (  IYl(l),  IY2(1),  NB  ) 

DO  M=2, LUMPS 

IF  (  NB  .EQ.  O  )  THEN 

CALL  DRAWB  (  1X2 (M) ,  IY2 (M) ,  1  ) 

ELSE 

CALL  DRAWB  (  IXl (M) ,  IYl (M) ,  1  ) 

END  IF 
END  DO 
ELSE 
END  IF 

DO  J=2, LUMPS 

DLIN  »  LIN  /  (LUMPS- 1) 

DL ( J)  =  (J-l)  *  DLIN 
END  DO 

IF  (  TTI  .GT.  0.0  )  THEN 
CALL  NEWTIME  (  TTI , 1  ) 

CALL  NEWTIME  (  TI(K),0  ) 

IF  (  NBED  .EQ.  NB  )  THEN 

CONTINUE 

ELSE 

CALL  FRAME2  (  NB,1  ) 

CALL  FRAME 2  (  NBED.O  ) 

END  IF 
ELSE 

CALL  NEWTIME  (  TI(K),0  ) 

CALL  FRAME2  (  NBED , O  ) 

END  IF 
TTI  =  TI(K) 

IF  (NOUTl  .EQ.  4)  GO  TO  120 
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WRITE  (lO,-)  PH  ,  NBED  , T , NT ,  NY  (K) 

WRITE (10,110)  TI(K) 

GO  TO  130 

120  WRITE (4 ,110)  TI(K) 

HO  FORMAT ( '  BED  COORDINATE ’, T30 ,* MOLE  FRACTION  IN  BED#1  ’  , 

lT48,’M0LE  FRACTION  IN  BED#2 ' , 3X , ' (TIME  T= ’ , F8 . 2 , ’ SEC . ) ’ , / ) 
130  CONTINUE 

DO  L=l, LUMPS 
XCA1 (L , K)  =  CA1 (L) 

XCA2 (L , K)  =  CA2(L) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
XCB1 (L,K)  =  CBl (L) 

XCB2 (L , K)  =  CB2 (L) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IF  (  NBED  . EQ .  O  )  THEN 
XAl(L.K)  =  XCA1(L,K)  /  Cl 
XA2 (L , K)  «  XCA3(L,K)  /  C2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
XBl (L , K)  «  XCBl (L , K)  /  Cl 
XB2(L,K)  =  XCB2 (L , K)  /  C2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

IXABSl (L , K)  =  276  *  IINC-(L-l) 

IXABS2 (L , K)  =  676  -  IINC-(L-l) 

IF  (  NOUTl  .EQ.  4  )  THEN 


(  XBl(L.K) 


IYABSl (L ,K)  ) 
IYABS2 (L , K)  ) 

IYABSl (L,K)  ) 
IYABS2(L,K)  ) 


CALL  SCALEl  (  XA1(L,K), 
CALL  SCALE2  (  XA2 (L , K) , 
ELSE 

CALL  SCALEl  (  XBl(L.K), 
CALL  SCALE2  (  XB2 (L , K) , 
END  IF 
ELSE 

M  =  LUMPS-*- 1  -L 

XAl(M.K)  =  XCA1 (L , K)  /  Cl 

XA2(M,K)  »  XCA2 (L , K)  /  C3 

XBl(U.K)  =  XCBl (L,K)  /  Cl 
XB2 (M , K)  =  XCB2 (L , K)  /  C2 


IXABS2 (L , K)  =  276  +  IINC«(L-1) 

IXABSl (L.K)  -  676  -  IINC*(L-1) 

IF  (  NOUTl  .EQ.  4  )  THEN 

CALL  SCALEl  (  XA1(M,K),  IYABSl (M,K)  ) 
CALL  SCALE2  (  XA2(M,K),  IYABS2(M,K)  ) 
ELSE 

CALL  SCALEl  (  XBl(M.K),  IYAES1(M,K)  ) 
CALL  SCALE2  (  XB2(M,K),  IYABS2(U,K)  ) 

END  IF 
END  IF 
END  DO 

CALL  INITIALl  (  IYABSl (1,K),  IYABS2(1,K), 

DO  1=2, LUMPS 

IF  (  NBED  .EQ.  O  )  THEN 

CALL  DRAWF  (  IXABSl (I, K),  IYABSl (I, K) ,  0  ) 
1X1(1)  =  IXABSl (I, K) 

IYl (I)  =  IYABSl (I , K) 

ELSE 

CALL  DRAWF  (  IXABS2(I,K),  IYABS2(I,K),  O  ) 
1X2(1)  =  IXABS2 (I , K) 

IY2 ( I )  =  IYABS2 (I , K) 

END  IF 
END  DO 

DO  1=1, LUMPS 


NBED  ) 


v‘v'v  vVvV,/W',,’'vV'  *v  »'-VkV'JN" -»■ 


IWCVCW  fWSBWSWW 


.  1  '  *  **» 


11  ^  I 


I  ij*  ^  ■  M  M:vL»i^««LM.a»I.^W  I 


IF  (NOUTl  .  EQ.  4)  GO  TO  148 
GO  TO  140 

148  WRITE (4 , 300)  DL (I) ,  XA1 (I , K) ,  XA2(I,K) 

GO  TO  180 

140  WRITE ( 1 0 , 300) DL ( I )  ,  XBl(I,K),  XB2(I,K) 

300  FORMAT  (  5X ,  F8 . 3 ,  14X,  FlO.8,  18X,  FlO.6  ) 

180  CONTINUE 
END  DO 

CALL  INITIAL2  (  IYABS1(1,K),  IYABS2(1,K),  NBED  ) 

DO  266  N=2, LUMPS 

IF  (  NBED  .Eq.  O  )  THEN 

CALL  DRAWB  (  IXABS2(N,K),  IYABS2(N,K),  O  ) 

1X2 (N)  =  IXABS2 (N , K) 

IY2 (N)  =  IYABS2 (N , K) 

ELSE 

CALL  DRAWB  (  IXABSl(N,K),  IYABS1(N,K),  O  ) 

IX1(N)  =  IXABSl (N , K) 

IY1(N)  =  IYABS1 (N , K) 

END  IF 
NB  =  NBED 
266  CONTINUE 

XBR022  =  XABR  /  1.06 
DY  =  XBR022 
DPY  =  DY  -  280.0 
IDY  =  IFIX  (  DPY  ) 

IABS  =  412  -  IDY 
IF  (  K  -GT.  1  )  THEN 
WRITE (0,600)  IA 

600  FORMAT  (  *  P [806 ,  ’ 13  1  ]  *  ) 

IF(IABS.LT.IA)  THEN 
ID  =  IA  -  IABS 
WRITE (0 ,610)  ID 

610  FORMAT  (  ’  W(R ,  NO)  V  []  V  [806  ,  -  ’  13  ’  ]  ’  ) 

ELSE 

ID  =  IABS  -  IA 
WRITE (0,620)  ID 

620  FORMAT  (  ’  W(R , Nl ) V [] V [806 , + ’ 13 ’ ] ’  ) 

END  IF 
ELSE 

WRITE (0,630)  IABS 

630  FORMAT  (  ’  P [806 , 412] V [ , ’ 13 ’ ] ’  ) 

END  IF 
WRITE (0,636) 

636  FORMAT  (  ’  W(NO) ’  ) 

C  WRITE (8, 200) DL (I) ,XA1 (I,K) ,XA2(I,K) 

C  200  FORMAT (8X ,F8.3,11X,E10.4, 7X , ElO . 4) 

range(l)  =  O. 
range (2)  =  18 . 
range (3)  =  O. 

range (4)  =  1 . 
nxaxia  =  20 

encode  (nxaxia , 268 , xaxia) 

268  format  (  ’BED  COORDINATE  (IN.)’  ) 
nyaxia  =  13 

encode  (nyaxia  ,  380  ,  jraxia) 

380  format  (  ’MOLE  FRACTION’  ) 
nxaxie60  =  11 

encode  (nxaxiaSO , 268 , xax ia60) 

268  format  (  ’TIME  (SEC.)’  ) 
ny  ax i aSO  =  0 


v;/v;vV /'a'  /v‘.a 
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366 


lOOO 


B2C  O 

0220 

eeee 

0460 

0400 

0300 


encode  (nyaxisSO , 366 , yaxisSO) 
format  (  ’LSO  (IN.)’  ) 

IA  =  IABS 
CONTINUE 

IF  ( ( (N0UT1 . EQ . 6) . AND .  (UOD (NT , INC)  . EQ . O) ) . OR . ( (NOUTl . EQ . 6) 

.  AND  .  (AMOD  (T ,  TCYC)  .  EQ  .  O .  )  )  )  THEN 
NBO  =  NSO  +  1 
DO  M=l, LUMPS 
IF  (  NBED  .EQ.  O  )  THEN 
ERR (M)  =  (CB1  (U) /Cl) - . lO 
ELSE 

L  =  LUMPS+l-M 

ERR (M)  =  (CB2(L)/C2)  -  .10 

END  IF 

IF  (  ABS (ERR (M) )  .LE.  0.015  )  THEN 
L60(N60)  =  DL(M) 

ELSE 
END  IF 
END  DO 

TSO (N60)  =  T 

WRITE (11,*)  TSO (NBO) , L60 (NSO) , NBED , LUMPS 
ELSE 
END  IF 

IF  (  NT  .EQ.  NY (6)  )  THEN 
RANGESO (1)  =  O. 

RANGE60(2)  =  O. 

RANGESO (3)  =  O. 

RANGESO (4)  =  O. 

CALL  ZETAPLT (11, 1 , NSO , RANGESO , TSO , LSO , NSO+2 .XAXISSO , 
NXAXISSO , , YAXISSO , NYAXIS60) 
call  setaplt (10,6 , lumps , range , DL ,XBl , LUMPS+2 , xaxie , 
nxaxie ,y axle , ny axis) 

call  aetaplt (20,6 , lumps , range , DL , XB2 , LUMPS+2 , xaxis , 
□xaxie , yaxis , nyaxie) 

CLOSE  (  UNIT=0 ,  STATUS= ’ SAVE ’  ) 

REWIND  O 

OPEN  (  UNIT=0 ,  FILE= * REGIS ’ ,  STATUS= ’ OLD ’  ) 

CALL  TXTERASE  (O) 

CALL  REGISTART 
CALL  PLTERASE 

READ  (  0,  0220,  END-0000  )  TEMP 
WRITE  (3,*)  TEMP 
TYPE  0220, TEMP 
FORMAT  (  A  ) 

GO  TO  0200 
CALL  REGISOUT 

CLOSE  (  UNIT=0,  STATUS=’ DELETE’  ) 

WRITE (6, 0460)  NULL 

FORMAT  (A,  ’  PRESS  "C"  TO  CONTINUE  i  ’  ) 

READ (6, 0300)  KEYIN 
FORMAT  (  A  ) 

IF  (  KEYIN  .EQ.  CONTINUE  .OR.  KEYIN  .EQ.  LOWC  )  THEN 
CALL  TXTERASE  (O) 

CALL  REGISTART 
CALL  PLTERASE 
CALL  REGISOUT 
go  to  17 
ELSE 

GO  TO  0400 
END  IF 
ELSE 

GO  TO  6 
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END  IF 

IF  (  MOD (NT, INC)  )  6,6,6 
IM  =  IM  ♦  1 

XBROa(IM)  =  XABR  /  1.06 
IF  (  NOUTl  .  EQ.  1  )  THEN 
WRITE ( 1 , 84)  XBR02(IM),  T 
FORMAT  (  8X ,  F10.6,  26X,  FlO.3) 

82  IF  (  ND  . EQ.  0  )  THEN 
CALL  TXTERASE  (  266  ) 

CALL  REGISTART 
CALL  PLTERASE 
CALL  FRAME 
TYPE  80 

80  FORMAT  (  * ! P [ 130 , 330] 1 ,  t  ) 

GO  TO  47 
ELSE 

CALL  STIME  (T,  IX,  TF ,  600.,  130  ) 

CALL  SCALEM  (  XBR02(IM),  IY ,  1.0,  1.0  ) 

TYPE  86,  IX,  IY 

86  FORMAT  (  ’ !V[’I3' , ’13*]  * ,  I  ) 

END  IF 

ELSE  IF  (  NOUTl  .EQ.  2  )  THEN 
IF  (  ND  .EQ.  O  )  THEN 
CALL  TXTERASE  (  266  ) 

CALL  REGISTART 
CALL  PLTERASE 
CALL  FRAME 
TYPE  86 

86  FORMAT  (  ’ !P [130, 330] ’ ,  f  ) 

GO  TO  47 

ELSE 

CALL  STIME  (  T,  DC,  TF,  600. ,  130  ) 

CALL  SCALEM  (  WV(ND.l),  IY,  0.06,  0.1  ) 

TYPE  87,  IX,  IY 

87  FORMAT  (  ’ !V[ '13' , ’13’] ’ ,  S  ) 

WRITE (2 , 68)  WV (ND , 1) ,  WV(ND,2),  T 

68  FORMAT  (  8X,  FlO.6,  llX,  FlO.6,  lOX,  PlO.3  ) 

END  IF 

ELSE  IF  (  NOUTl  .EQ.  3  )  THEN 
IF  (  ND  .EQ  .0  )  THEN 
CALL  TXTERASE  (  266  ) 

CALL  REGISTART 
CALL  PLTERASE 
CALL  FRAME 
TYPE  88 

88  FORMAT  (  ’ !P[130,330] ’ ,  S  ) 

GO  TO  47 

ELSE 

CALL  STIME  (  T,  IX,  TF ,  600.,  1LJ  ) 

CALL  SCALEM  (  WV(ND,2),  IY,  0.08,  0.1  ) 

TYPE  89 ,  IX,  IY 

80  FORMAT  (  ’ !V[’I3’ , ’13’] ’ ,  I  ) 

WRITE (2 , 69)  WV (ND , 1 ) ,  WV(ND,2),  T 
60  FORMAT  (  8X ,  FlO.8,  llX,  FlO.6,  10X,  F10.6,  6X ,  FlO.3  ) 
END  IF 
ELSE 
GO  TO  82 
END  IF 
47  CONTINUE 
ND=ND+ 1 

IF  (  PH- 180 .  )  2,3,3 
WV  (ND ,  1 )  =W  (3) 


TSTEP2  .AND .  T  .LE.  TP  )  THEN 


WV(ND,2)=W(1) 

GO  TO  4 
WV(ND,1)=-W(l) 

WV(ND,2)=-W(3) 

XV (ND , 1) =XBR02 (IR) 

XV  (ND  ,  2)  =0  . 

XV  (ND  , 3) =1 . 

TV  (ND)  =T 
NT=NT+ 1 
T  =  NT  »  DT 

IF  (  T  .LT.  TSTEPl  .AND.  T  .La.  TF  )  THEN 
WBRLOLD  =  WBRL 
GO  TO  1 

ELSE  IF  (  T  .GE.  TSTEPl  .AND.  T  .LT.  TRTEP2  )  THEN 
WBRL  =  BRSTEP 
GO  TO  1 

ELSE  IF  (  T  .GE.  TSTEP2  -AND.  T  .LE.  TP  )  THEN 
WBRL  =  WBRLOLD 
GO  TO  1 
ELSE 

WBRL  =  WBRLOLD 
GO  TO  18 
END  IF 

CONTINUE 
CALL  REGISOUT 
WRITE (6, 42)  NULL 

FORUAT  (A,  1  PRESS  "C"  TO  CONTINUE  :  ’  ) 

READ (6 ,41)  KEYIN 
FORUAT  (  A  ) 

IF  (  KEYIN  Eq.  CONTINUE  .OR.  KEYIN  .EQ.  LOWC  )  THEN 
CALL  TXTERASE  (O) 

CALL  REG I START 
GO  TO  43 
ELSE 

GO  TO  44 
END  IF 

CALL  PLTERASE 
CALL  REGISOUT 
range2(l)  =  O. 
rtng«2 (2)  =  O. 

rtng«2(3)  =  O. 
range2  (4)  *  1  . 
r»ng«3(l)  =  O. 
rtng<3(3)  »  O. 
range3(3)  *  -.Ol 
rtng<3(4)  =  .04 
nxtxiiS  =  ii 

encode  (  nxxxi*2  ,  257,  xaxia2  ) 

format  (  ’TIRE  (SEC.)’  ) 
nyaxia2  =  16 

•ncode  (  nyaxiaS,  368,  yaxia2  ) 
format  (  ’02  MOLE  FRACTION’  ) 
nxaxiaS  =  11 

•ncode  (  nxaxiaS,  862,  xaxia3  ) 

format  (  ’TIRE  (SEC.)’  ) 

nyaxia3  =  18 

encode  (  nyaxiaS,  963,  yaxia3  ) 

format  (  'BREATHING  FLOWRATE’  ) 

call  aetaplt  (  30,  3,  nd ,  range2,  tv,  xv ,  nd+2,  xaxia2, 
nxaxiaS,  yaxia2,  nyaxiaS  ) 

call  eetaplt  (  40,  2,  nd ,  ranges,  tv,  wv ,  nd  +  2,  xaxia3, 
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nxkxia3,  ]raxi*3 ,  nyaxia3  ) 


lOO  FORMAT  (  14X,  8E12.3  )  , 

GO  TO  17  i 

COO  WRITE (6,010)  I 

OlO  FORMAT (/////// . . ’,/////,  1 

1 ’  PROGRAM  TERMINATED  BY  OPERATOR ’ , 

2 //////• . ’,//////)  ; 

END 


BLOCK  DATA 

REAL  LIN , KA , KB , KC 

COMMON/SYST/DT , PSUP , POUT , TF , WBRL , VS , VO , IIN , TCYC 
COMMON/GEOM/DBYIN , DVSIN , DVOIN , LIN , DIAOIN , DIAIIN 
COMMON/CAP/KA , KB , KC ,BB , BC , D 
COMMON/VOID/E 

DATA  DT,  PSUP,  POUT,  TP,  WBRL,  VS,  VO,  IIN,  TCYC  /.Ol, 

&  40.0,  16.0,  30.0,  10.0,  0.0,  0.0,  O,  10.7  / 

DATA  DBYIN,  DVSIN,  DVOIN,  LIN,  DIAOIN,  DIAIIN  / 
a  .076,  .306174,  .43866,  16.6,  6.73,  2.18  / 

DATA  KA,  KB,  KC,  BB,  BC ,  D  /  .2133,  .007607,  .06223, 

&  .6414,  .6266,  200.  / 

DATA  E  /  .37  / 

END 

SUBROUTINE  MENU 

REAL  LIN ,KA,KB , KC ,LINl , KA1 ,KBl ,KC1 

DIMENSION  TI (6) ,  NY<8) ,  TI4(6) ,  NY4<6) ,  TI6(6),  NY6 (6) 

BYTE  KEYIN 
BYTE  YES 
BYTE  loyea 
BYTE  NO 
BYTE  lono 

COMMON/SYST/DT , PSUP , POUT , TF , WBRL , VS , VO , IIN , TCYC 
COMMON/GEOM/DBYIN , DVSIN , DVOIN , LIN , DIAOIN , DIAIIN 
COMMON/BEDV/LUMPS , DZ , L , D1 , D2 , MA , MB , MC , RT , AI , AO , AMWOl , 

AMWL 1 , AMW02 , AMWL2 , WABR , WBBR 
COMMON/VOID/E 

COMMON/CAP/KA , KB , KC , BB , BC , D 
COMMON/NFINISH/NEXIT 

COMMON/STEP/NSTEP ,TSTEP1 , TSTEP2 , BRSTEP 
C0MM0N/NPRINT/N0UT1 ,TI,NY 
DATA  YES, loyea/88, 121/ 

DATA  NO , lono/78 , 1 10/ 

NEXIT=0 
NSTEP=0 
lO  CALL  SYSTl 
CALL  GEOUl 
CALL  CAP! 

WRITE (7, 107) 

WRITE (6, 107) 

107  FORMAT C  ENTER  CORRESPONDING  PARAMETER  #  TO  CHANGE  PARAMETER’,/, 
1’  ENTER  "O"  TO  SIGNIFY  PARAMETERS  ARE  CORRECT- -RUN  SIMULATION’,/ 
2’  ENTER  "60"  TO  INSERT  A  STEP  CHANGE  IN  BREATHING  FLOWRATE’,/, 

4’  ENTER  "OB"  TO  EXIT  PROGRAM’) 

N  =  O 
READ  • , N 

IF  (  N  .  EQ.  0  )  GOTO  600 
IF  (  N  .EQ.  6B  )  GOTO  600 
IF  (  N  . EQ.  80  )  GOTO  OOO 
EPSI  =  .OOOOOOl 
WRITE (7, 120) 

WRITE (6, 120) 

120  FORMAT  (///  ,  ’ . . . . 
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1  //,’  ENTER  THE  PARAMETER  AND  PRESS  "RETURN"’,/, 

2//T40, ’CURRENT  VALUE’,/) 

IP  (  N  .LT.  8  )  GOTO  lOO 
IF  (  N  .LT.  12  )  GOTO  200 
IF  (  N  .LE.  16  )  GOTO  300 
GOTO  800 
lOO  CONTINUE 

IF  (  N-2  )  169,156,156 
166  IF  (  N-3  )  164,171,171 

171  IF  (  N-4  )  170,175,174 

169  WRITE (8, 160)  PSUP 
WRITE (7, 160)  PSUP 

180  FORMAT (’  SUPPLY  PRESSURE  (PSIA) ’ .T40.F8.4, /) 

READ  -,PSUP1 

IF (PSUP1 . LT . ESPI)  GOTO  lO 

PSUP=PSUP1 

GOTO  lO 

184  WRITE (8 , 168)  POUT 

188  FORMAT ( ’  OUTLET  PRESSURE  (PSIA) 1 ,T40,F8.4,/) 

READ  • .POUTl 

IF (POUTl . LT . EPSI)  GOTO  10 

POUT=POUTl 

GOTO  lO 

170  WRITE (8, 172)  TF 

172  FORMAT ( ’  FINAL  OBSERVATION  TIME  (SEC) ’ ,T40 , F8 . 4 , /) 

READ  * , TFl 

IF(TFl.LT.EPSI)  GOTO  lO 

TF=TFl 

GOTO  lO 

176  WRITE (8, 177)  TCYC 

177  FORMAT ( ’  CYCLE  TIME  (SEC) ’ , T40.F8 . 4 , /) 

READ  • ,TCYCl 

IF (TCYC1 . LT . EPSI)  GOTO  lO 

TCYC=TCYCl 

GOTO  lO 

174  WRITE (8, 178)  WBRL 

178  FORMAT ( ’  BREATHING  MASS  FLOWRATE  (STD  LIT/MIN) ’ ,T40,F8.4, /) 
READ  - , WBRLl 

IF ( WBRL 1 .LT. EPSI)  GOTO  178 
WBRL= WBRLl 
GOTO  lO 
178  GO  TO  lO 
200  CONTINUE 

IF (N-7)  209,216,206 

206  IF (N-9)  226,236,207 

207  IF(N-IO)  236,246,266 

209  WRITE (8, 210)  DBYIN 

210  FORMAT ( ’  BY-PASS  VALVE  DIAMETER  (IN) ’ , T40 , F8 . 4 , /) 

READ  * , DBYIN 1 

IF (DBYINl . LT . EPSI)  GOTO  lO 

DBYIN=DBYIN1 

GO  TO  lO 

216  WRITE (8, 220)  DVSIN 

220  FORMAT (’  SUPPLY  VALVE  DIAMETER  (IN) ’ , T40 , F8 . 4 , /) 

READ  * , DVSINl 

IF (DVSIN1 . LT . EPSI)  GOTO  10 

DVSIN=DVSINl 

GOTO  lO 

226  WRITE (6, 230)  DVOIN 

230  FORMAT ('  OUTLET  VALVE  DIAMETER  (IN) ’ , T40, F8 . 4 , /) 

READ  * , DV0IN1 

IF (D VO INI . LT . EPSI)  GOTO  lO 
DVOIN=DVOINl 


GOTO  lO 

336  WHITE (6 , 340)  LIN 

340  FORMAT (’  BED  LENGTH  (IN) ’ , T40 , F8 . 4 , /) 

READ  * ,LINl 

IF (LINl . LT . EPSI)  GOTO  lO 

LIN=LIN1 

GOTO  lO 

346  WRITE (6, 360)  DIAOIN 

360  FORMAT (’  OUTER  BED  DIAMETER  (IN) ’ , T40 , F8 . 4 , /) 

READ  * , DIA0IN1 
IF (DIA0IN1.LT. EPSI)  GOTO  lO 
DIA0IN=DIA0IN1 
GOTO  lO 

366  WRITE (6, 360)  DIAIIN 

360  FORMAT (’  INNER  BED  DIAMETER  (IN) ’ ,T40 , F8 . 4 , /) 

READ  • , DIAIIN1 
IF (DIAIIN1.LT. EPSI)  GOTO  lO 
DIAIIN=DIAIIN1 
GOTO  lO 
300  CONTINUE 

IF (N- 13) 300 , 316 , 306 
306  IF (N- 16) 336, 336, 346 
300  WRITE (6, 310)  KA 
310  FORMAT  (  *  KA’,  T40,  F8.4,  /  ) 

READ  - ,  KA1 

IF  (  KA1  .LT.  EPSI  )  GOTO  lO 
KA  =  KA1 
GOTO  lO 

316  WRITE (6, 330)  KB 

330  FORMAT  (  *  KB ’ ,  T40,  F8 . 4 ,  /  ) 

READ  - ,  KB1 

IF  (  KB1  .LT.  EPSI  )  GOTO  lO 
KB  =  KB1 
GOTO  lO 

336  WRITE (6, 330)  KC 

330  FORMAT  (  ’  KC ’ ,  T40 ,  F8 . 4 ,  /  ) 

READ  *,  KC1 

IF (  KC1  .LT.  EPSI  )  GOTO  lO 
KC  =  KCl 
GOTO  lO 

336  WRITE (8, 340)  D 

340  FORMAT ( ’  DIFFUSION  COEFFICIENT  (> lOO . ) ’ , T40 , F8 . 4 , /) 

READ  •,D1 

IF (Dl . LT . EPSI)  GOTO  lO 

D=Dl 

GOTO  lO 

346  WRITE (6, 360)  E 

360  FORMAT (’  VOID  FRACTION  « 1 . ) ' , T40 , F8 . 4 , /) 

READ  • , El 

IF (El .LT. EPSI)  GOTO  lO 

E=El 

GOTO  lO 

600  WRITE (8, 610)  TF 
WRITE (7 , 610)  TF 

610  FORMAT (////  ’ . . . . '//. 

1 ’  ENTER  THE  TIME  AT  WHICH  THE  STEP  CHANGE ’ , / 

3’  IN  BREATHING  FLOWRATE  IS  TO  OCCUR’,//, 

3'  THE  FINAL  OBSERVATION  TIME  IS  CURRENTLY  ’,F7.3,’  SECONDS’,///) 
NSTEP=1 
READ  « , TSTEPl 
616  WRITE (6 ,618)  TF 
WRITE (7,618)  TF 

618  FORMAT  (////  ’ . . . *//. 
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1*  ENTER  THE  TIME  AT  WHICH  THE  STEP  CHANGE’,/ 

3'  IN  BREATHING  FLOWRATE  IS  TO  END’,//, 

3’  THE  FINAL  OBSERVATION  TIME  IS  CURRENTLY  ’,F7.2,’  SECONDS’,///) 
READ  * , TSTEP2 

WRITE (7, 620)  TSTEP 1 , TSTEP2 , WBRL 
WRITE (6, 620)  TSTEP 1 , TSTEP2 , WBRL 

620  FORMAT  (//  ’ . . . . . ’//  , 

1 ’  ENTER  THE  NEW  BREATHING  FLOWRATE  OCCURING  AS ’ , / , 

2*  A  STEP  CHANGE  FROM  T=  ',F7.2,’  TO  T=  ’,F7.2,’  SECONDS’,//, 

3’  THE  CURRENT  BREATHING  FLOWRATE  IS’,F8.2,*  STD  LIT/MIN ’ , ///) 

READ  • , BRSTEP 

WRITE (7, 630)  BRSTEP, TSTEPl , TSTEP2 
WRITE (8, 630)  BRSTEP, TSTEP 1 .TSTEP2 
630  FORMAT (////, 

1*  . ’,////. 

2’  THE  BREATHING  FLOWRATE  OF  ’,F8.2,’  STD  LIT/MIN’,/ 

3’  WILL  BE  INPUT  AS  A  STEP  AT  T=  ' ,F7 . 2,  *  SECONDS’/, 

3’  AND  WILL  END  AT  T=  ’,F7.2,  ’  SECONDS’////, 

4’  . ’////) 

GOTO  lO 
800  CONTINUE 

WRITE(7 ,810) 

WRITE (8,810) 

810  FORMAT  (/ ///// ,  . *,/////, 

. ’  THE  FOLLOWING  PLOTS  ARE  AVAILABLE  TO  YOU  ON  THE  TERMINAL’,///, 

.*  1  MOLE  FRACTION  OF  OXYGEN  VS. TIME  (DATA  FILE=MOLEFRAC . DAT) ’ , / , 
.’  2  INLET  MASS  FLOWRATE  VS.  TIME  (DATA  FILE=FLOWRaTE . DAT) ’ , / , 

.*  3  OUTLET  MASS  FLOWRATE  VS.  TIME  (DATA  FILE=FLOWRATE . DAT)  ’  ,  /  , 

.'4  A  DYNAMIC  SIMULATION  FOR  ONE  CYCLE  (OXYGEN  VS.  TIME,’,/, 

DATA  FILE=PROFILE . DAT) * ,/, 

. ’  6  DYNAMIC  SIMULATION  OF  PROPAGATION  OF  CONTAMINANT  ALONG ’ , / , 
THE  BED  LENGTH ’ , / , 

•///. 

. ’  PRESS  THE  CORRESPONDING  #  TO  HAVE  THE  OUTPUT  SHOWN  AS  THE ’ , / , 

.’  SIMULATION  TRANSPIRES’ ,/////) 

READ (6 , 830) NOUTl 
830  FORMAT (II) 

IF  (NOUTl  .  EQ.  4  )  GOTO  760 
IF  (NOUTl  . EQ.  6  )  GOTO  780 
GO  TO  400 
760  DO  J=1 ,6 

FJ  =  FLOAT (J) 

TI4(J)  =  2. -FJ  -  1 . 

NY4 ( J)  =  IFIX  (  TI4 ( J) /DT  ) 

END  DO 

TI4 (8)  =  TCYC 

NY4 (8)  =  IFIX  (  TI4 (8) /DT  ) 

DO  M=1 ,8 
TI(M)  =  TI4 (M) 

NY (M)  =  NY4(M) 

END  DO 
GO  TO  400 
80  DO  J=2 , 6 

FJ  =  FLOAT ( J) 

TI6 ( J)  =  (  TF/6.  )  *  (FJ-1.) 

NY 6 ( J )  =  IFIX  (TI6(J)/DT  ) 

END  DO 

TI6(1)  =  1.0 

NY6 ( 1 )  =  IFIX  (  TI6(1)/DT  ) 

TI6 (8)  =  TF 

NY 6(8)  =  IFIX  (  TI6 (8) /DT  ) 

DO  M= 1 , 8 
T I (M)  =  TI6 (M) 


END  DO 

400  WRITE (6,410) 

WRITE (7 , 410) 

410  FORUAT (////, ’  YOU  HAVE  THE  FOLLOWING  OPTIONS  : 

./////, TlO, *1 .  ENTER  "O"  TO  RUN  SIMULATION.’,/, 

. TlO , '  2  .  ENTER  "19"  TO  CHANGE  THE  PARAMETERS.’,/, 

. TlO , ’ 3 .  ENTER  "99"  TO  EXIT  PROGRAM.’,/////, 

ENTER  THE  CORRESPONDING  #  AND  PRESS  "RETURN"’,/////) 

READ  - ,  N69 

IF(N8W.EQ.O)  GOTO  830 
IF (N89 . EQ . 19)  GOTO  lO 
IF (NS9 . EQ . 99)  GOTO  900 
800  CONTINUE 

WRITE (7 ,810) 

WRITE (8 ,810) 

810  FORMAT  (/////  ,  ' . *,/////, 

1’  THE  PARAMETER  VALUE  WAS  NOT  INPUT  CORRECTLY’,/ 

3’  TRY  AGAIN’ ,/////,  ’ . ',////) 

GOTO  lO 

820  WRITE (6, 826) 

WRITE (7, 836) 

826  FORMAT (////, ’  DO  YOU  WISH  TO  DOUBLE  THE  NUMBER  OF  SPACE  LUMPS  IN’/, 
1’  THE  SIMULATION  FOR  BETTER  ACCURACY  ?  (  CURRENT  NUMBER  USED  ’/, 

3’  IS  60  )  (Y/N) ’,/////) 

READ  (6,830)  KEYIN 
830  FORMAT  (  A  ) 

IF  (  (KEYIN  .E(J.  YES)  .OR.  (KEYIN  .EQ  ■  loye»)  )  THEN 

LUMPS  =  lOl 

ELSE 

LUMPS  =61 
END  IF 
GO  TO  lOOO 
900  CONTINUE 
NEXIT  =  1 
lOOO  RETURN 
END 

SUBROUTINE  SYST1 

COMMON/SYST/DT , PSUP , POUT , TF , WBRL , VS , VO , IIN , TCYC 
WRITE (8, 360)  PSUP , POUT , TF , TCYC , WBRL 
WRITE (7, 360)  PSUP , POUT , TF , TCYC , WBRL 
360  FORMAT (’  CURRENT  SYSTEM  PARAMETERS  ARE:’,/, 

1’  1 . ’ , T6 , ’  SUPPLY  PRESSURE-’ ,T30,F8.3,T40, ’PSIA’ ,/, 

3’  3 .  ’ , T6 ,  ’  OUTLET  PRESSURE-’ ,T30,F8.2,T40,  ’PSIA’ ,/, 

3’  3 . ’ , T6 , '  FINAL  OBSERVATION  TIME- ’, T30 , F8 . 3 , T40 ,’ SEC ’,/ , 

4'  4 .  ’ , T6 ,  ’  CYCLE  TIME- ’ , T30 , F8 . 2 , T40 ,  ’ SEC ’ , / , 

6'  6 . ' , T6 , ’  BREATHING  FLOWRATE- ', T30 , F8 . 2 , T40 ,’ STD  LIT/MIN ’ ) 

RETURN 

END 

SUBROUTINE  GE0M1 
REAL  LIN 

COMMON/GEOM/DBYIN , DVSIN , DVOIN ,LIN,DIAOIN,DIAIIN 
WRITE (8 ,460)  DBYIN, DVSIN, DVOIN, LIN, DIAOIN , DIAIIN 
WRITE(7 ,460)  DBYIN, DVSIN, DVOIN, LIN, DIAOIN, DIAIIN 
460  FORMAT (’  CURRENT  GEOMETRIC  PARAMETERS  ARE:’,/, 

1’  8 . ’ , T6 , '  BY-PASS  VALVE  DIAMETER  =  ’ , T30 ,F8 . 4 , T40 , ’ IN ’ , / , 

2’  7 . ’ , T6 , ’  SUPPLY  VALVE  DIAMETER  =  ’ , T30 , F8. 4 , T40 , ’ IN ’ , / , 

3’  8 . ’ , T6 , ’  OUTLET  VALVE  DIAMETER  =  ’ , T30 , F8 . 4 , T40 , ’ IN ’ , / , 

4’  9 . ' , T6 , ’  BED  LENGTH  =  ’ ,T30 , F8 . 4 , T40 , 'IN’,/, 

6’  lO . ’ , T6 , ’  OUTER  BED  DIAMETER  =  ’ ,T30 , F8 . 4 , T40 , ’ IN ’ , / , 

8’  11.’, T6 , ’  INNER  BED  DIAMETER  *  ' , T30 , F8 . 4 , T40 , ’IN’) 


01  *  W  10 


RETURN 

END 

SUBROUTINE  CAPl 
REAL  KA , KB , KC 

COMMON/CAP/KA ,  KB ,  KC ,  BB  ,  BC ,  D 

COMMON/VOID/E 

WRITE (7 , 660) KA , KB , KC , D ,E 

WRITE (6 , 660) KA , KB , KC , D , E 

660  FORMAT (’  THE  CURRENT  BED  PARAMETERS  AREt’,/, 


1  ’ 

12.  ’  , 

,  T6 ,  ’ 

KA  = 

’ , T30 , F8 . 4 , T40 , * KGMOLES 

02 

ADSORBED/KGMOLE 

GAS ’ , / , 

2  1 

13.  ’  , 

,  T6 ,  • 

KB  = 

' , T30 , F8 . 4 , T40 , ’ KGMOLES 

CO 

ADSORBED /KGMOLE 

GAS ’ , / , 

3  ' 

14.  ’  , 

i  T6 ,  ’ 

KC  = 

’ , T30 , F8 . 4 , T40 , ’ KGMOLES 

N2 

ADSORBED/KGMOLE 

gas • , / , 

4’ 

16.  ’ 

.  T6 ,  ’ 

DIFFUSION  COEFFICIENT  =  ' , T30 , F8 . 2 , T40 , ’ 1/SEC ’ , / 

» 

6  ’ 

16  .  ’  . 

,  T6 ,  * 

VOID 

FRACTION  =  ’ , T30 , F8 . 3) 

RETURN 

END 

SUBROUTINE  TIMEF (CVA , PH , PS , T) 

C ..... . CALCULATES  CVA(T)  AND  PS (T) 

DIMENSION  CVA (3), PS (2) 

DATA  PHI , PHI ,PH2 , PH3 , PH4 , PH6 , PH6/10 . ,20. ,160. ,180. , 
1  200 . , 340 . , 360 . / 

COMMON/SYST/DT , PSUP , POUT , TF , WBRL , VS , VO , IIN , TCYC 
COMMON/HILOBED/NBED , NB 

C . * . * . 

c 

C  PH  =  PHASE  ANGLE  OF  CYCLE  (DEGREES) 

C  TCYC  =  CYCLE  TIME  (SEC) 

C 

. . * . . . . . 

1  =  0 

PH  =  AMOD  (  360..T/TCYC  +  PHI,  360.  ) 

IF  (  PH-PHl  )  1,2,2 
1  RA  =  PH  /  PHI 

GO  TO  11 

IF  (  PH-PH2  )  3,4,4 
RA  =  1  . 

GO  TO  11 

IF  (  PH-PH3  )  6,6,6 
RA  =  (PH3-PH)  /  (PH3-PH2) 

GO  TO  11 

6  IF  (  PH-PH4  )  7,8,8 

7  RA  =  (PH-PH3)  /  (PH4-PH3) 

I  =  1 

GO  TO  11 

8  IF  (  PH-PH6  )  0,10,10 

0  RA  =  1  . 

I  =  1 
GO  TO  11 

lO  RA  =  (PH8-PH)  /  (PH8-PH6) 

I  =  1 

1  1  NBED  =  I 

IF  (  I  )  12,12,13 

12  PSfl)  =  PSUP 

PS (2)  =  POUT 
CVA ( 1 )  =  VS-RA 
CV  A ( 3 )  =  VQ.RA 

RETURN 

13  PS(1)  =  POUT 
PS (2)  =  PSUP 
CVA(l)  =  VO.RA 
CVa (3)  =  VS-RA 


-*.\V  ,v 


>  /  / 
'  vV. 


•  Hi  **•  • 
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RETURN 

END 

SUBROUTINE  INITIAL  (Cl,  C2 ,  Nl,  N2 ,  XA ,  XB,  P  ) 

REAL  L , MA , MB , MC , NAOl , NA02 , KA , KB , KC , Nl , N2 , NBOl , NB02 , LIN , NCOl , NC02 
DIMENSION  Nl (300,3) ,N2(300,3) ,XA(3) ,XB(3) ,P(2) 

COMMON/BEDV/LUMPS , DZ , L , D1 , D2 , MA , MB , MC , RT , AI , AO , AMWOl , AMWL1 , AMW02 
1  , AMWL2 , WBR 

COMMON/CAP/KA , KB , KC , BB , BC , D 
COMMON/VOID/E 

COMMON/GEOM/ DBYIN , DVSIN , DVOIN , LIN , DIAOIN , DIAIIN 


RT  =  GAS  CONSTANT. TEMP  (  N-M  /  KGMOLE  ) 

AITOT , AOTOT  =  INLET  AND  OUTLET  CONCENTRIC  BED  AREAS  (  M-«2  ) 
MA ,  MB  =  KGMOLECULAR  WEIGHTS  (  KG  /  KGMOLE  ) 

O.L  =  SUBSCRIPTS  FOR  Z=0  AND  Z=L 

AMW  =  AVERAGE  MOLECULAR  WEIGHT  (  KG  /  KGMOLE  ) 

1,2  =  SUBSCRIPTS  FOR  BED  1  AND  BED  2 
E  =  VOID  FRACTION 


RT  =  8314.  .  300. 

AITOT  =  3 . 14150  •  (DIAIIN*. 0254)*. 2  /  4. 

AOTOT  =  3.14160  -  (DIAGIN*»2-DIAIIN-.2)  -.0264.-2  /  4 
L  =  .0264  -  LIN 

PRINT  . , ’  BED  LENGTH  (M)  =  ’ ,  L 
PRINT  INLET  BED  AREA  (M.-2)  =  *,  AITOT 

PRINT  .,’  OUTLET  BED  AREA  (M--2)  =  *,  AOTOT 
PRINT  - , ’  VOID  FRACTION  =  * ,  E 
MA  =  32. 

MB  =  28. 

MC  =  28. 

AI  =  AITOT. E 
AO  =  AOTOT-E 

AMWOl  =  MA.XA(l)  +  MB-XB(l)  +  MC*  (1  .  -XA  (1 )  -XB (1)  ) 
AMWLl  =  MA-XA (2)  +  MB*XB(2)  +  MC* (1 . -XA(2) -XB(2) ) 
AMW02  =  AMWLl 

AMWL2  =  MA-XA (3)  +  MB«XB(3)  +  MC* (1 . -XA(3) -XB (3) ) 

Cl  =  P (1 ) *8 . 806E3  /  RT 
C2  =  P (2) *0 . 805E3  /  RT 


Cl -XA ( 1) 

C2-XA (3) 

Cl-XB(l) 

C2-XB(3) 

Cl-CAl-CBl 
C2-CA2-CB2 
CA1  /  KA 
CA2  /  KA 

(CBl/KB)  /  (1 . +  CB1  / (KB-BB) ) 
(CB2/KB)  /  (1 . +CB2/ (KB-BB) ) 
(CCl/KC)  /  (1 . -CC1/(KC-BC) ) 
(CC2/KC)  /  (1 . +  CC2/ (KC-BC) ) 


CA1  =  Cl -XA ( 1) 

CA2  =  C2-XA (3) 

CBl  =  Cl-XB(l) 

CB2  =  C2-XB (3) 

CCl  =  Cl-CAl-CBl 
CC2  =  C2-CA2-CB2 
NAOl  =  CA1  /  KA 
NA02  =  CA2  /  KA 
NBOl  =  (CBl/KB)  /  (: 
NB02  =  (CB2/KB)  /  ( 
NCOl  =  (CCl/KC)  /  ( 
NC02  =  (CC2/KC)  /  ( 
DO  1  1=1, LUMPS 
Nl  (I , 1)  =  NAOl  +  NB< 
Nl (I , 2)  =  NAOl 
Nl (1,3)  =  NBOl 
N2 (1,1)  =  NA02  ♦  NB< 
N2 (1,2)  =  NA02 
N2 (1,3)  =  NB02 
DZ  =  L/FLOAT (LUMPS- 1 ) 
RETURN 
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SUBROUTINE  BEDS  (  P,  W,  XA,  XB,  XABR,  XBBR ,  CVA ,  PS,  DTG , 
CAX  ,  CA2,  CBl,  CB2 ,  Cl ,  C2  ) 


-SOLVES  BED  EQUATIONS 


REAL  L , Nl , N2 , UA , MB , MC 

DIMENSION  CAX (300) , CA2 (300) , CBl (300) , CB2 (300) , Nl (300 , 3) , N2 (300 , 3) , 
X  Cl (300,3) , C2 (300 , 3) ,XAT(2) ,XBT(2) ,XA0(2) ,XB0(2) 

1  , W (3) , XA (3) , XB (3) ,CVA(3) ,PS(2) ,PN(2) ,P(2) 

COMMON/BEDV/LUMPS , DZ , L , Dl , D2 , UA , MB ,MC , RT , AI , AO , AMWOl , AMWLl , AMW02 
1  , AMWL2 , WBR 

COMMON/VOID/E 

COMMON/SYST/DT , PSUP , POUT , TF , WBRL , VS , VO , IIN , TCYC 

COMMON/CAP/KA , KB , KC , BB , BC , D 

COUMON/TIME/T 


C  =  TOTAL  GAS  STREAM  CONCENTRATION  (  KGMOLE/  M--3  ) 
CA  =  GAS  A  STREAM  CONCENTRATION  (  KGMOLE  /  M--3  ) 

N  =  ADSORBED  PRASE  CONCENTRATION  (  KGMOLE  /  M--3  ) 
G  =  GAS  PHASE  CONCENTRATION  (  KGMOLE  /  M--3  ) 

U  =  STREAM  VELOCITY  (  M  /  SEC  ) 

Y  =  MOLECULAR  FLUX  (  KGMOLE  /  U- -2-SEC  ) 

L  =  BED  LENGTH  (  M  ) 

D  =  DIFFUSION  COEFFICIENT  (  1  /  SEC  ) 

BETA  =  D-AREA  RATIO  (  1  /  SEC  ) 

DT  =  TIME  SrEP  (  SEC  ) 

LUMPS  =  NUMBER  OF  SPACE  LUMPS 

1,2  =  SUBSCRIPTS  FOR  BED  1  AND  BED  2 

VBR  =  MIXING  VOLUME  FOR  BREATHING  FLOW  (  M--3  ) 

BR  =  SUBSCRIPT  INDICATING  BREATHING  VARIABLE 

CP  =  PLENUM  MOLE  CONCENTRATION  (  KGMOLE  /  U--3  ) 

VP  =  BREATHING  PLENUM  VOLUME  (  M--3  ) 


DT=DTG 
DO  4  1=1 , 2 

PN (I)  =  PS(I) -8.806E3 
IF (IIN . EQ . 1) GO  TO  1 
IIN  =  1 

CALL  INITIAL  (Cl,  C2,  Nl ,  N2 ,  XA,  XB ,  P  ) 

VPMUL  =  1. 

VP  =  .0018  *  VPMUL 

PRINT  PLENUM  VOLUME  (M--3)  =  ’ , VP 

DO  1=1,7 
PRINT  - 
END  DO 

WRITE (8, lO) NULL 

FORMAT  (  A,  ’  PLEASE  WAIT . *  ) 

VBR  =  .6- VP  +  .0003 

C?  =  . 96- AMAXl (  PN ( 1 ) , PN (2)  )  /  RT 

DO  3  1=1 , 2 

XAT(I)  =  XABR 

XBT(I)  =  XBBR 

Dl  =  D 

D2  =  D 

WBR  =  WBRL-4 . 72E-5  /  2.206 

CALL  CONCEN  (Cl , C2 , CA1 , CA2 , CBl , CB2 , W, Gl , G2 , XA , XB , Nl , N2 , PN , CVA) 
XABR  =  XAT  (2) 

XBBR  =  XBT (2) 

AM WBR  =  MA-XABR  +  MB -XBBR  +  MC- ( 1 . -XABR -XBBR) 


DO  6  1=1,2 
XAO(I)  =  XAT (I) 
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v\' 


XBO(I)  =  XBT(I) 
DO  2  1=1, LUMPS 


Nl (I , 1 )  =  Nl(I.l) 
Nl (I , 2)  =  Nl (I , 2) 
Nl (1,3)  =  Nl (1,3) 
N2 (1,1)  =  N2 (1,1) 
N2 (1,2)  =  N2 (1,2) 
N2 (1,3)  =  N2 (1,3) 
CONTINUE 
FAC  =  CP  *  VBR  * 
TAUI  =  WBR  /  FAC 
XAT(l)  =  XAO(l)  + 
XBT(l)  =  XBO(l)  + 
XAT (2)  =  XAO (2)  + 
XBT (2)  =  XBO (2)  + 
P(l)  =  Cl-RT  /  8. 
P (2)  =  C2-RT  /  8. 
AMWOl  =  Ka-XA(I) 
AMWL1  MA-XA(2) 
AMW02  =  AMWLl 
AMWL2  =  MA-XA(3) 
RETURN 
END 


Dl*(Gl(I,l)-Cl) *DT 
Dl-  (Gl (1,2) -CA1 (I)  )  -DT 
Dl- (Gl (1,3) -CBl  (I) ) -DT 
D2*  (G2 (1,1) -C2) -DT 
D2- (G2 (I , 2) -CA2 (I) ) -DT 
D2- (G2 (1 , 3) -CB2  (I) ) -DT 


AMWBR 

TAUI- (XA (2) -XAO (1 )  )  -  DT 
TAUI-  (308(2)  -XBO  (1)) -DT 
TAUI- (XAO (1 ) -XAO (2) ) -DT 
TAUI- (XBO (1) -XBO (2)) -DT 
89SE3 
886E3 

+  MB-XB(l)  +  MC*  (1 .  -XA(1)  -XB  (1)  ) 
f  MB-XB (2)  +  MC- (1 . -XA(2) -XB (2) ) 

+  MB-XB  (3)  +  MC-  ( 1 .  -XA  (3)  -XB  (3)  ) 


SUBROUTINE  CONCEN (CGI , CG2 , CA1 , CA2 , CBl , CB2 , W, Gl , G2 , 

XA , XB , N 1 , N2 , PN , CVA) 

SOLVES  FOR  CONCENTRATIONS 
REAL  L , N 1 , N2 , UA , MB , MC 

DIMENSION  CA1 (300) , CA2 (300) ,CB1 (300) , CB2 (300) ,Nl (300,3) ,N2(300,3) 
, Gl (300,3) , G2 (300 , 3) 

DIMENSION  Gil (300) ,GI2(300) ,PN(2) ,W(3) ,XA(3) ,XB(3) ,CVA(3) 
C0MM0N/BEDV/LUMPS,DZ,L,Dl,D2,MA, MB, MC,RT,AI, AO, AMWOl, AMWLl, AMW02 
, AMWL2 , WBR 
COMMON/VOID/E 

C0MM0N/VBEQ/W1 ,W2,W3,U01 ,ULl ,U02,UL2,PNl ,PN2,VA1 ,VA2,VA3 
, BGI 1L , BGI1L2 , BGI2L , BGI2L2 , BL21 , BL22 , PI ,P2 


BED  PRESSURE  (  N  /  M--2  ) 

SUPPLY  PRESSURES  (  N  /  M--2  ) 

STREAM  VELOCITY  (  M  /  SEC  ) 

MASS  FLOWRATE  (  KG  /  SEC  ) 

DISCHARGE  COEFFICIENT  TIMES  VALVE  AREA  (  M--2  ) 


TOTAL  GAS  STREAM  CONCENTRATION 
GAS  "A"  STREAM  CONCENTRATION 
GAS  "B"  STREAM  CONCENTRATION 


(  KGMOLE  /  M--3  ) 
(  KGMOLE  /  M* » 3  ) 

(  KGMOLE  /  M--3  ) 


GAS  "C"  STREAM  CONCENTRATION  (  1 
BED  DISTANCE  COORDINATE  (  M  ) 
MOLECULAR  FLUX  OF  GAS  "A"  =  CA-U 
MOLECULAR  FLUX  OF  GAS  "B"  =  CB-U 
MOLE  FRACTION  OF  GAS  A 
MOLE  FRACTION  OF  GAS  B 


VA1  =  CVA ( 1 ) 

VA2  =  CVA (2) 

VA3  =  CVA (3) 

BETA1  =  (l.-E)-Dl  /  E 
BETA2  =  ( 1 . -E) -D2  /  E 
PN 1  =  PN ( 1 ) 

PN2  =  PN(2) 

BL21  =  BETA1 - . 6-L 
BL22  =  BETA2- . 6-L 


(  KGMOLE  /  M' 


(  KGMOLE  /  M* -2-SEC  ) 
(  KGMOLE  /  M- -2-SEC  ) 


vvvvv  V 
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Cl  =  CGI 
C2  =  CG2 

--CALCULATE  G  AND  INTEGRAL  OF  G 
Gil  (1)  =0.0 
GI2 ( 1 )  =  0.0 
DO  1  1=1, LUMPS 
CALL  CAPZ  (  Nl,  Gl,  I  ) 

CALL  CAPZ  (  N2,  G2 ,  I  ) 

DO  2  1=2, LUMPS 
IU1  =  1-1 

GIl(I)  =  GIl(IMl)  -  DZ-(  Gl(I,l)+Gl (IMl , 1 )  )  /  2. 

GI2 (I)  =  GI2 (IMl )  *  DZ-  (  G2 (I , 1) +G2 (IMl , 1)  )  /  2. 

LUM2  =  (LUMPS* 1)  /  2 
BCIlLa  =  BETA1  -  Gil (LUM2) 

BGI1L  =  BETA1-GI1  (LUMPS)  -  BGHL2 
BGI2L2  =  BETA2  -  GI2 (LUM2) 

BGI2L  =  BETA2-GI2 (LUMPS)  -  BGI2L2 
--SOLVE  BED/VALVE  EQUATION  FOR  C 
CALL  VAQ  (Cl,  C2  ) 

W(1)=W1 
W(2)=W2 
W(3) =W3 
CG1=C1 
CG2=C2 

--SOLVE  FOR  CONCENTRATION  OF  GAS  A  AND  B 
CALL  CONCAB  (  Cl , CA1 , CBl , Gl , Gil , UOl , ULl , XA (1)  , XA (2)  , 

1  XB(1) ,XB(2) ,BETA1 ,AI/AO  ) 

CALL  CONCAB  (  C2 , CA2 , CB2 , G2 , GI2 , U02 , UL2 , XA (2) , XA (3) , 

1  XB(2) ,XB(3) ,BETA2,A0/AI  ) 

RETURN 

END 

SUBROUTINE  CONCAB  (  CG ,  CA,  CB,  G,  Gl ,  UOG ,  ULG,  XAO,  XAL , 

1  XBO,  XBL,  BETAG,  ARG  ) 

.--SOLVES  FOR  CONCENTRATION  OF  GAS  A  AND  B 
REAL  L2,L 

DIMENSION  CA (300) ,CB(300) ,G(300,3) ,GI(300) ,XA(2) ,XB(2) , 

1  U(300) ,YA(300) ,YB(300) 

C0MM0N/BEDV/LUMPS,DZ,L,D1 ,D2,MA,MB,MC,RT,AI,A0,AMW0l ,AMWL1 ,AMW02 
1  , AMWL2 , WBR 
COMUON/VOID/E 
COMMON/TIME/T 
C  =  CG 

BETA  =  BETAG 
UL  =  ULG 
UO  =  UOG 
XA ( 1 )  =  XAO 
XA  (2)  =  XAL 
XB ( 1 )  =  XBO 
XB  (2)  =  XBL 
AR  =  ARG 

LUM2  =  (LUMPS* 1)  /  2 
GIL2  =  GI(LUM2) 

L2  =  . 6-L 

IF  (  UO* -2  +UL--2  )  22,22,23 


■ZERO  FLOW  SOLUTION 


DO  24  1=1, LUMPS 
CAd)  =  G  (1 , 2) 
CB ( I )  =  G (I  ,  3) 
RETURN 
CONTINUE 


/. Vv'AV.  A  .•  /. 
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CA ( 1 )  =  C-XA(l) 

CB ( 1 )  =  C-XB(l) 

U(l)  =  UO 

YA (1)  =  UO-CA(l) 

YB  (1 )  =  UO-CB(l) 

CA (LUMPS)  =  C-XA(2) 

CB (LUMPS)  =  C-XB (2) 

YA (LUMPS)  =  UL-CA (LUMPS) 

YB (LUMPS)  =  UL-CB (LUMPS) 

Z  =  0.0 

-ICR  IS  INDEX  INDICATING  U=0  IN  BED 
ICR  =  1 

ICRl  =  ICR  -  1 

-SOLVE  FOR  U  IN  BED 

DO  3  1=2 .LUMPS 
IMl  =  1-1 
Z  =  2  +  DZ 

IF  (  IM1-LUM2  )  36,34,34 

U(I)  =  U(l)  -BETA* (  GI(I)/C  -  Z  ) 

GO  TO  36 

U(I)  =  AR-U (LUM2)  +  BETA- (  (GI (I) -GIL2) /C  -  Z  +  L2  ) 
CONTINUE 

-CHECK  SIGN  CHANGE  OF  U 

IF  (  U (I) -U (IMl )  )  8,8,3 
ICR  =  I 
ICRl  =  ICR-1 

-CALCULATE  INTERIOR  INITIAL  FLUX  (  U0<0 ,  UL>0  ) 


ULAST  =  U (IMl ) 

IF  (  IMl  ,E(J.  LUM2  )  ULAST=AR. ULAST 
DZCR  =  DZ-U(I)  /  (  U(I) -ULAST  ) 

DTC  =  DZCR  -  DZ 

Rl  =  U (ICR)  /  DZCR 

R3A  =  (  G  (I,  3) -G  (IMl  ,2)  )  /  DZ 

R3B  =  (  G (I, 3) -G (IMl, 3)  )  /  DZ 

R2A  =  G (IMl , 2)  -  R3A-DTC 

R2B  =  G (IMl , 3)  -  R3B-DTC 

Al A  =  BETA-R2A-R1  /  (  BETA+Rl  ) 

AlB  =  BETA-R2B-R1  /  (  BETA+Rl  ) 

A2A  =  BETA-R3A-R1  /  (  BETA+2.-R1  ) 

A2B  =  BETA-R3B-R1  /  (  BETA+2.-R1  ) 

YA (ICR)  =  A1A-DZCR  +  A2A -DZCR -DZCR 
YB (ICR)  =  A1B-DZCR  +  A2B- DZCR -DZCR 
YA (ICRl )  =  A1A-DTC  +  A2A-DTC-DTC 
YB  (ICRl )  =  AlB -DTC  +  A2B-DTC-DTC 
CA ( ICR)  =  YA (ICR)  /  U (ICR) 

CB (ICR)  =  YB (ICR)  /  U (ICR) 


CA (ICRl)  =  YA (IC 
CB (ICRl )  =  YB (IC 
GO  TO  3 
ICR  =  I 
ICRl  =  I 
YA (ICR)  =  0.0 
YB  (ICR)  =  0.0 
IF  (  UO  )  3,26,3 


YA (ICRl )  /  ULAST 
YB (ICRl )  /  ULAST 


f  >*  JV.  (V 
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ICR1  =  1 
YA ( 1 )  =0.0 
YB (1 )  =  0.0 
CONTINUE 
UL  =  U  (LUMPS) 

••LOGIC  FOR  FORWARD  INTEGRATION 
ISIG  =  1 

IF  (  UL  )  14,14,15 
1ST  =  1 

IF  (  UO  )  12,20,21 

CONTINUE 

IT  =  ICR-2 

IF  (  IT  )  12,12,11 

ICR  =  1 

ICRl  =  1 

GO  TO  12 

CONTINUE 

1ST  =  ICR 

IT  =  LUMPS-ICR 

IF  (  IT  )  12,12,11 

CONTINUE 


•SOLVE  FOR  CA 


DO  lO  J=1 , IT 
I  =  1ST  +  ISIG* J 
IM1  =  I-ISIG 

IF  (  IM1-LUM2  )  33,30,33 
IF  (  ISIG  )  31,33,32 
YA(IMl)  =  YA(IMl)  /  AR 
YB(IUl)  =  YB(IUl)  /  AR 
U(IM1)  =  U(IM1)  /  AR 
GO  TO  33 

U(IU1)  =  AR.U(IMl) 

YA(IMl)  =  AR»YA(IM1) 

YB(IMl)  =  AR-YB(IMl) 

CONTINUE 
RO  =  U(IMl) 

IF  (  I  . EQ.  LUM2  .AND.  ISIG  .Eft.  -1  )  U(I)=AR*U(I) 

R1  =  (  U (I) -RO  )  /  DZ 
R2A  =  G (IMl , 2) 

R2B  =  G (IMl , 3) 

R3A  =  (  G (I , 2) -R2A  )  /  DZ 
R3B  =  (  G  (1 , 3) -R2B  )  /  DZ 
IF  (  Rl  )  8,7,8 

AOA  =  (  R2A-R3A*RO/BETA  )  *  RO 
AOB  =  (  R2B-R3B*RO/BETA  )  *  RO 
A1A  =  R3A*R0 

AlB  =  R3B.RO 

BR  =  BETA  /  RO 

YA (I)  =  (  YA (IMl ) -AOA  ) -EXP (-BR*DZ)  +  AOA  +  A1A.DZ 

YB (I)  =  (  YB (IMl ) -AOB  ) *EXP (-BR-DZ)  +  AOB  +  A1B*DZ 

GO  TO  lO 

A2A  =  BETA-R3A-R1  /  (  BETA+2.*Rl  ) 

A2B  =  BETA.R3B.R1  /  (  BETA+2.-R1  ) 

A1A  =  (  BETA- (R3A-RO+R2A-R1  )  -  2.-A2A-RO  )  /  (  BETA+Rl  ) 

AlB  =  (  BETA- (R3B-RO+R2B-R1  )  -  2.-A2B-RO  )  /  (  BETA+Rl  ) 

AOA  =  (  BETA-RO-R2A-A1A-RO  )  /  BETA 

AOB  =  (  BE' A-RO-R2B-A1B-RO  )  /  BETA 

IF  (  RO  )  26,27,20 

YA (I)  =  A1A-DZ  +  A2A-DZ-DZ 

YB (I)  =  AlB-DZ  +  A2B-DZ-DZ 

GO  TO  28 


A1A-DZ 

AlB-DZ 
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CONTINUE 
BR  =  BETA/R1 

YA ( I )  =  (  YA (IMl ) -AOA  )-(  l.+Rl-DZ/RO  )--(-BR) 
+  A2A-DZ-DZ 

YB (I)  =  (  YB  (IMl ) -AOB  )-(  1 . +R1 -DZ/RO  )--(-BR) 
+  A2B-DZ-DZ 

CONTINUE 

IF  (  U(I)  )  18,19,18 

CA  (I)  =  C»G  (1 , 2)  /  G  (1 , 1 ) 

CB ( I )  =  C-C (1,3)  /  G (I , 1) 

GO  TO  lO 

CA  ( I )  =  YA (I )  /  U(I) 

CB ( I )  =  YB (I)  /  U(I) 

CONTINUE 

IF  (  ICR  .ESI.  1  )  CO  TO  13 
IF  (  DZ  )  13,13,13 
CONTINUE 


A1 a-dz 


A  IB*  DZ 


■logic  for  backward  integration 


DZ  =  -DZ 
ISIG  =  -1 

IF  (  UL  )  10,18,17 

1ST  =  LUMPS 

IT  =  LUMPS  -  ICR 

IF  (  ICR  Eq.  ICRl  )  IT=IT- 1 

IF  (  IT  )  13,13,11 

CONTINUE 

1ST  =  ICRl 

IT  =  ICRl - 1 

IF  (  IT  )  13,13,11 

CONTINUE 

XA ( 1 )  =  CA ( 1 )  /  C 
XA(3)  =  CA (LUMPS)  /  C 
XB(1)  =  CB ( 1 )  /  C 
XB(2)  =  CB (LUMPS)  /  C 
DZ  =  FLOAT (ISIG)  -  DZ 
XAO  =  XA ( 1 ) 

XAL  =  XA  (2) 

XBO  =  XB  (  1 ) 

XBL  =  XB  (2) 

RETURN 

END 

SUBROUTINE  CAPZ  (  N,  G,  I  ) 

REAL  KA  ,  KB  ,  KC  ,  N  ,  NA  ,  NB  ,  NC 
DIMENSION  N ( 300 , 3 ) , G ( 300 , 3 ) 

COMMON/CAP/KA , KB , KC , BB , BC , D 
NA  =  N(I , 2) 

NB  =  N  (I  ,  3) 

NC  =  N  (1 , 1 )  -  NA  -  NB 

G ( I , 2)  =  KA-NA 

G  (I , 3)  =  KB-NB  /  (l.-NB/BB  ) 

GC  =  KC*NC  /  (l.-NC/BC  ) 

G  (1 , 1 )  =  G  (1 , 2)  ♦  G  (1 , 3)  »  GC 

RETURN 

END 

SUBROUTINE  VALVE  (  A ,  PU ,  PD ,  W  ) 

** -CALCULATES  VALVE  MASS  FLOWRATE  KG/SEC 
REAL  KI 

DATA  KI,  Cl,  C2 ,  CK  /  .714,  .00900,  .00234,  280  / 

PUL  =  PU 

PDL  =  PD 

PR  =  PDL  /  PUL 
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SIG  =  1 . 

IF  (  1 . -PR  )  4,3,2 
W  =  O  . 

RETURN 
4  PR  =  1 .  /  PR 

PUL  =  PDL 
SIG  =  -1  . 

3  CONTINUE 

IF  (  PR  .LT.  .628  )  GO  TO  1 
PRR  =  1.  -  PR--CK 

W  =  SIG-A-Cl-PUL- (PR-.KI) -SQRT(PRR) 

RETURN 

1  W  =  SIC* A-C2-PUL 

RETURN 
END 

SUBROUTINE  VAQ  (  CGI ,  CG2  ) 

C-. ....ITERATION  AND  FALSE  POSITION  TO  SOLVE  VALVE/BED  EQUATIONS 
EXTERNAL  Fl ,F3 

COkUON/BEDV/LUUPS , DZ , L , Dl , D2 , UA , MB , MC , RT , AI , AO , AMWOl ,  AlfWLl , AMW02 
1  , AMWL2.WBR 

COUMON/VBEQ/Wl , W2 , W3 , UOl , UL1 , U02 , UL2 , PNl , PN2 , VA1 , VA2 , VA3 
1  , BGI 1L , BGIIL2 , BGI2L , BGI2L2 , BL21 ,BL22,P1 ,P2 
Cl  =  CGI 
C2  =  CG2 
Pi  =  Cl «RT 
POl  =  Pi 

J  =  1 

to  FC2  =  F2  (  C2  ) 

1  CN2  =  C2- (  1.-FC2  ) 

IF  (  ABS (CN2-C2) /C2  .Le.  l.E-8  )  GO  TO  2 
FCN2  =  F2  (  CN3  ) 

IF  (  FC2-FCN2  )  4,3,3 

3  C3  =  CN2 
FC2  =  FCN2 
GO  TO  1 

4  C2  -  FP  (  F2 ,  FCN3,  CN2 ,  FC2 ,  C2  ) 

2  CONTINUE 

FCl  =  Fl  (  Cl  ) 

6  CN 1  =  Cl  * (  1 . -FCl  ) 

IF  (  ABS (CNl -Cl ) /Cl  .LE.  l.E-8  )  GO  TO  8 
FCNl  =  Fl  (  CNl  ) 

IF  (  FCl* FCNl  )  8,7,7 

7  Cl  =  CNl 
FCl  =  FCNl 
GO  TO  6 

8  Cl  =  FP  (Fl,  FCNl,  CNl,  FCl,  Cl  ) 

8  IF  (  ABS (Pl-P01)/Pl  .LE.  l.E-8  )  GO  TO  8 

J  =  J»1 
POl  =  Pi 

IF  (  J  LE  20  )  GO  TO  lO 
PRINT  lOO 

loo  FORMAT  (  IX,  ’.FAILED  TO  CONVERGE-’  ) 

8  CGI  =  Cl 

CG2  =  C2 
RETURN 
END 

FUNCTION  Fl  (  CGI  ) 

C . BED  l  VALVE/BED  EQUATION 

COMMON/BEDV/LUMPS , DZ , L , Dl , D2 , UA , MB , MC , RT , AI , AO , Atfwoi , AMWLl , AMW02 
1  ,  A1TWL2  ,  WBR 

COMMON'/ VBEQ/Wl  ,W2,W3,U01  ,  UL  1  ,  U02  ,  UL  2  ,  PN 1  ,  PN2  ,  VA1  ,  VA2  ,  VA3 
1  . BG 1 1 L , BGI 1 L2 . BGI2L , BGI2L2 , BL2 1  , BL22.Pl ,P2 

Cl  =  CGI 
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PI  =  Cl -RT 

CALL  VALVE  (  VA1 ,  PN1 ,  PI,  Wl  ) 

CALL  VALVE  (  VA2 ,  Pi,  P2,  W2  ) 

IF  (  P1-P2  )  1,2,3 
WBRl  =  O. 

GO  TO  4 
WBRl  =  . 6-WBR 
GO  TO  4 
WBRl  =  WBR 
CONTINUE 

UOl  =  Wl  /  (  Cl-AMWOl-AI  ) 

ULl  =  (  W2+WBR1  )  /  (  Cl-AMWLl-AO  ) 

AR  =  AI  /  AO 

FI  =  1.  -  (  AR-BGHL2+BGI1L  )  /  (  (UL1+BL21-AR- (U01-BL21) ) -Cl  ) 

RETURN 

END 

FUNCTION  F2  (  CG2  ) 

-.-BED  2  VALVE/BED  EQUATION 

COMMON /BEDV /LUMPS , DZ , L , D1 , D2 , MA , MB , MC , RT , AI , AO , AMWOl , AMWL1 , AMW02 
1  , AMWL2 , WBR 

COMMON /VBEQ/Wl , W2 , W3 , UOl , ULl , U02 ,  UL2 , PNl ,PN2, VA1 , VA2.VA3 
1  , BGI1L , BGX1L2 , BGI2L , BGI2L2 , BL21 , BL22 , Pi ,P2 
C2  =  CG2 
P2  =  C2-RT 

CALL  VALVE  (  VA2 ,  PI,  P2,  W2  ) 

CALL  VALVE  (  VA3 ,  P2 ,  PN2 ,  W3  ) 

IF  (  P2-P1  )  1,2,3 
WBR2  =  O. 

GO  TO  4 
WBR2  =  .6-WBR 
GO  TO  4 
WBR2  =  WBR 
CONTINUE 

U02  *  <  W2-WBR2  )  /  (  C2-AMWO2-A0  ) 

UL2  =  W3  /  (  C2-AMWL2-AI  ) 

AR  =  AO  /  AI 

F2  =  1.  -  (  AR-BGI2L2+BGI2L  )  /  (  (UL2+BL22-AR- (U02-BL22) ) -C2  ) 

RETURN 

END 

FUNCTION  FP  (  F,  FXR,  XR,  FXL ,  XL  ) 

--USES  FALSE  POSITION  TO  FIND  ZERO 
DATA  EPS  /  l.E-fl  / 

1=0 

X  =  XL  -FXL* (XL-XR)  /  (  FXL -FXR  ) 

FX  =  F  (  X  ) 

I  =  1-1 

IF  (  ABS (FX)  .LT.  EPS  )  GOTO  4 
IF  (  I  .GT.  60  )  GOTO  4 
IF  (  FX-FXL  )  1,1,2 
XR  =  X 

FXL  =  FXL -FXR  /  (  FXR+FX  ) 

FXR  =  FX 
GOTO  3 
XL  =  X 

FXR  =  FXR -FXL  /  (  FXL+FX  ) 

FXL  =  FX 
GOTO  3 
FP  =  X 
RETURN 
END 

subroutine  zetiplt (nplot , ncva , npte , range , x , y , ny , xax is , nxaxis , 

yaxis , nyaxis) 


dimension  x (300) , y (302 , 1 1) , c (8) , range (4) ,xaxie(10) ,yaxis(10) 

data  axlenx , axleny  /  7.0,  5.0  / 

call  plots  (  O,  0,  nplot  ) 

call  plot  (  0.5,  0.6,  -3  ) 

call  factor  (  0.66  ) 

if  (  (range (1) . eq.O. )  .and.  (range (2) . eq. O . )  )  then 
call  scale  (  x,  axlenx,  npta,  1  ) 
else 

x(npts+l)  =  range(l) 

x(npts+2)  =  (range (2) -range (1) ) /axlenx 
end  if 

if  (  (ra,.^e(3)  .eq.O.)  .and.  (range  (4)  .  eq .  0  .  )  )  then 
call  scale  (  y,  axleny,  npts ,  1  ) 
else 

y(npts+l,l)  =  range (3) 

y(npte+2,l)  =  (range (4) -range (3) ) /axleny 
end  if 

if  (  ncvs  .ge.  2  )  then 
do  n  =  2, ncvs 

y (npts+1 , n)  =  y(npts+l,l) 
y(npts+2,n)  =  y(npts+2,l) 
end  do 
end  if 

call  axis  (  0.0,  0.0,  xaxis,  -nxaxis,  axlenx,  0.0, 
x(npts+l),  x(npts+2)  ) 

call  axis  (  0.0,  0.0,  yaxis,  nyaxis,  axleny,  00.0, 
y(npts+l,l),  y (npts+2, 1)  ) 
do  m  =  1 , ncvs 

call  line  (  x,  y(l,m),  npts,  1,  O,  2  ) 
end  do 

call  plot  (  0.0,  0.0,  000  ) 

return 

end 

SUBROUTINE  SCALEl  (  Y,  YABS  ) 

COMMON/HILOBED/NBED 
INTEGER  YABS 
DY  =  1 . O  -  Y 
DPY  =  DY  *  160.0 
IDY  =  IF IX  (  DPY  ) 

YABS  =  IDY  +  lOO 

RETURN 

END 

SUBROUTINE  SCALE2  (  Y,  YABS  ) 

COMMON/HILOBED/NBED 
INTEGER  YABS 
DY  =  1.0  -  Y 
DPY  =  DY  *  1600 
IDY  =  IF IX  (  DPY  ) 

YABS  =  IDY  +  300 

RETURN 

END 

SUBROUTINE  INITIAL1  (  YABS1 ,  YABS 2 ,  NBED  ) 

INTEGER  YABSl , YABS2 
IF  (  NBED  .EQ.  O  )  THEN 
WRITE (0 , lOO)  YABSl 
ELSE 

WRITE (0,100)  YABS2 
END  IF 

lOO  FORMAT  (3X ,  ’  PfaTS.’IS’]’  ) 

RETURN 

END 

SUBROUTINE  INITIAL2  (  YABSl,  YABS2 ,  NBED  ) 
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INTEGER  YABSl , YABS2 
IF  (  NBED  .EQ.  O  )  THEN 
WRITE (9,100)  YABS2 
ELSE 

WRITE (9,100)  YABSl 
END  IF 

lOO  FORMAT  (3X,  ’  P[576,’I3']'  ) 

RETURN 

END 

SUBROUTINE  DRAWF  (  XABS,  YABS ,  NM  ) 

C  DRAW  VECTORS  FROM  PREVIOUS  CURSOR  POSITIONS  TO  CURRENT  POINTS 
C  FORWARD  PLOTTING 

INTEGER  XABS,  YABS 
WRITE (9,100)  NM, XABS, YABS 
lOO  FORMAT  (3X,  *  W (R , N ’ II ’ ) V  [] V  [  • 13 1  ,  *  13  ’  ]  *  ) 

RETURN 

END 

SUBROUTINE  DRAWB  (  XABS,  YABS,  NM  ) 

C  BACKWARD  PLOTTING 

INTEGER  XABS, YABS 
WRITE (9,100)  NM,  XABS,  YABS 
lOO  FORMAT  (3X ,  ’  W(R , N ’ II ’ ) V [] V [ • 13 ’ ,  ’13*]’  ) 

RETURN 

END 

SUBROUTINE  TXTERASE  (NTICKS) 

C  ERASE  THE  VT126  TEXT  SCREEN 
LOGICAL  ERASE 
DIMENSION  ERASE (4) 

DATA  ERASE  /  27 ,  ’ [’ ,  1 2 ’ ,  ’ J ’  / 

IF  (  NTICKS. NE.O  )  THEN 
DO  1=1,7 

CALL  DELAY  (NTICKS) 

END  DO 
ELSE 
END  IF 

TYPE  lO,  ERASE 
lO  FORMAT  (  ’ + 1 ,  6A1 ,  8  ) 

RETURN 

END 

SUBROUTINE  REGISTART 
C  SETS  VT126  INTO  GRAPHICS (REGIS)  MODE 
LOGICAL  GRAPHIC 
DIMENSION  GRAPHIC (3) 

DATA  GRAPHIC  /  27 ,  ’P*,  ’p’  / 

TYPE  lO,  GRAPHIC 
lO  FORMAT  (  ’+ ’ ,  6A1 ,  8  ) 

RETURN 

END 

SUBROUTINE  PLTERASE 
C  ERASES  PLOTS 

LOGICAL  ERASE 
DIMENSION  ERASE (4) 

DATA  ERASE  /  'S’,  ’(’,  ’E’,  ’)’  / 

TYPE  lO,  ERASE 
lO  FORMAT  (  •+’,  6A1,  8  ) 

RETURN 

END 

SUBROUTINE  FRAMEl 
WRITE (9, llO) 

llO  FORMAT  (  •  P [250 , 6] T (S2) "OBOGS  SIMULATION"  ) 

WRITE(9, 112) 

112  FORMAT  (  1  P [30, 30] T (Si) "TIME  =  SEC.  ) 

WRITE(9, 114) 
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114  FORMAT  (  ’  P[260,06]T"l"P[aSOIa4fi]T"0,,  ’  ) 

WRITE (8 , 118) 

118  FORMAT  (  *  P [280 , 306] T" 1 "P [280 , 445] T"0* ’  ) 

WRITE (8 , 120) 

120  FORMAT  (  ’  P[400,80]T*BED  #1"’  ) 

WRITE(9 , 130) 

130  FORMAT  (  *  P [ 400 , 280 ] T " BED  #2"'  ) 

WRITE (9, 140) 

140  FORMAT  (  ’  P [30 , 256] T"INLET" ’  ) 

WRITE (9 , 160) 

160  FORMAT  (  *  P [6 , 276] T" PRESSURE" ’  ) 

WRITE (9, 160) 

180  FORMAT  (  *  P [ 196 , 1 12] T*R . V . ■ *  ) 

WRITE (9, 170) 

170  FORMAT  (  ’  P [96 , 430] T"EXHAUST" ’  ) 

WRITE (9 , 180) 

180  FORMAT  (  ’  P  [116,460]T*GAS" '  ) 

WRITE (9 , 1 90) 

190  FORMAT  (  ’  P [824 , 288] T"P . 0 . » ’  ) 

WRITE (9 , 192) 

192  FORMAT  (  *  P [880 , 1 12] 1 "B . P . "  ’  ) 

WRITE (9 , 193) 

193  FORMAT  (  ’  P [720 , 412] VOOOOOP [730 , 403] T"0*» ’  ) 

WRITE (9 , 196) 

196  FORMAT  (  ’  P [720 , 368] VOOOOOP [730 . 349] T* 30** ’  ) 

WRITE (9 , 198) 

198  FORMAT  (  ’  P [720, 132] VOOOOOP [730, 123] T" 10O%" ’  ) 

WRITE (9, 300) 

200  FORMAT  (  ’  P [276, lOO ]V [676, 100]V [676, 260] V [376, 260] V [276,100]  •) 

WRITE (9 , 210) 

210  FORMAT  (  ’  P [276 , 300] V [876 , 300] V [676 , 460] V [276 , 460] V [276 , 300] *  ) 

WRITE  (9, 300) 

300  FORMAT  (  ’  P [676 , 162 ] V [ 670 , 182] V [670 , 132 ] V [720 , 1 32] V [720 ,412] '  ) 

WR I TE  ( 9 , 310) 

310  FORMAT  (  ■  V [ 670 ,412]V[670, 382 ] V [676 , 382] P [ 67 6 , 1 82 ]  ’  ) 

WRITE (9, 400) 

400  FORMAT  (  ’  V [600, 182] V [600, 267] V [606 ,372] V [600, 277] V [600,363] *  ) 

WRITE (9,410) 

410  FORMAT  (  '  V [ 67 6 , 362] P [620 , 1 82 ] V [670 ,182]V[670, 362 ] V [620 , 362]  '  ) 

WRITE (9, 460) 

460  FORMAT  (  ’  V [ 620 , 277 ] V [6 1 6 , 272] V [620 , 267] V [62C , 182 ] ’ ) 

WRITE ( 9 , 460) 

460  FORMAT  (  ’  P [720 , 282] V [760 , 262] V [760 , 282] V [720, 282] ’  ) 

WRITE (9 , 600) 

600  FORMAT  (  ’  P [276 , 162] V[240, 162] V[240, 133] V[180, 132] ’  ) 

WRITE ( 9 ,610) 

610  FORMAT  (  *  V [ 1 80 , 1 62[ V [ 1 1 8 , 1 63] V [ 1 1 8 , 362] '  ) 

WRITE (9 , 620) 

620  FORMAT  (  ’  P [ 27  6 , 1 82] V [ 240 , 1 82 ] P [ 1 80 , 182]V[138, 182]V[138, 382 ] 

1  '  ) 

WRITE ( 9 , 640) 

640  FORMAT  (  '  P[11B,282]V[118, 420] V [ 1 38 , 420] V[138,382J  '  ) 

WRITE ( 9 , 600) 

800  FORMAT  (  *  V [ 180 , 382 J  V [ 1 80 , 41 2 ] V [240 , 41 3] V [240 , 382 ] V [ 276 , 382]  ’  ) 

WRITE (9 ,810) 

810  FORMAT  (  '  P [ 1 80 , 282 J V [80 , 262] V [80 , 282] ’  ) 

WRITE ( 9 , 700) 

700  FORMAT  (  ’  V [ 180 , 383] V [ 180 , 362] V [ 138 , 362] V [ 138 , 282] P [ 276 ,362]  '  ) 

WRITE (9,710) 

710  FORMAT  (  ’  P [276 ,362] V [340, 362] V [240 , 182]P [180, 182] V [ 180, 262] ’  ) 
WRITE ( 9 , 730) 

730  FORMAT  (  ’  P [ 830 , 1 63 ] V [630 , 187 ] V [640 , 1 63]  ’  ) 

WR ITE (9 , 736) 
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736  FORMAT  (  ’  P [030 , 183] V [630 , 177] V [640, 182] ’  ) 

WRITE (8,740) 

740  FORMAT  (  ’  P [630 , 362] V [630 , 367] V [640 , 362] ’  ) 

WRITE (8, 760) 

760  FORMAT  (  ’  P [630 , 382] V [630 , 377] V [640 , 382] ’  ) 

RETURN 

END 

SUBROUTINE  NEWTIME  (  TIME , NW  ) 

WRITE (8,100)  NW.TIME 

lOO  FORMAT  (  ’  W(N’Il ’ ) P [1 10 , 30] T (Si ) • ’F6. 2 ’ ■ ’  ) 

RETURN 

END 

SUBROUTINE  FRAME2  (  NBED ,NW  ) 

IF  (  NBED  .EQ.  1  )  THEN 
WRITE (8,710) NW 

710  FORMAT  (  ’  W(N’Il ’ ) P [180 , 262] V [240 , 382] ’  ) 

ELSE 

WRITE (8 , 720)  NW 

720  FORMAT  (  ’  W (N 1 II ’ ) P [180 , 262] V [240 , 162] ’  ) 

END  IF 

WRITE (8, 730) NW 

730  FORMAT  (  '  W (N ’ 11 ’ ) P [630 , 182] V [630 , 167] V [640 , 1 62] ’  ) 

WRITE (0, 736) 

736  FORMAT  (  *  P [830 , 182] V [630 , 177] V [640, 182] ’  ) 

WRITE (8, 740) NW 

740  FORMAT  (  ’  W (N ’ II ’ ) P [630 , 362] V [630 , 367] V [640 , 362] '  ) 

WRITE (8 , 760) 

760  FORMAT  (  ’  P [830 , 382] V [830 , 377] V [640 , 382] ’  ) 

IF  (  NBED  . EQ.  1  )  THEN 
WRITE (8 , 800) NW 

800  FORMAT  (  ’  W(N’I1 ’ ) P [180 , 282] V [240 , 382]  ’  ) 

WRITE (8, 806) 

806  FORMAT  (  1  P [240 , 182] V [180 , 162] P [240, 182] V [180 , 182] ’  ) 
WRITE (8 , 810) NW 

810  FORMAT  (  ’  W (N ’ 11 ’ ) P [113 , 272] V [143 , 272] ’  ) 

WRITE (8,816) 

816  FORMAT  (  ’  P [136 , 288] V [] V [143 , 272] V [138 , 278] ’  ) 

WRITE (8, 820) NW 

820  FORMAT  (  *  W (N ’ 11 ’ ) P [680 , 372] V [626 , 372] ’  ) 

WRITE (8, 826) 

826  FORMAT  (  ’  P [820 , 388] V [] V [826 , 372] V [820 , 376] ’  ) 

WRITE (8, 830) NW 

830  FORMAT  (  ’  W (N ’ II ’ )P [610 , 320] V [610 , 280] ’  ) 

WRITE (8, 836) 

836  FORMAT  (  ’  P [607 , 286] V[]V[610, 280] V[613,286] ’  ) 

WRITE (8 , 840) NW 

840  FORMAT  (  ’  W(N’Il ’ )P [810, 172] V [686 , 172] ’  ) 

WRITE (8, 846) 

846  FORMAT  (  ’  P [680 , 168] V [] V [686 , 172] V [680 , 176] ’  ) 

WRITE (8, 860) NW 

860  FORMAT  (  ’  W(N ’ 11 ’ ) P [ 128 , 362] V [ 128 , 382] ’  ) 

WRITE (8, 866) 

866  FORMAT  (  ’  P [ 124 , 386] V [] V [128 , 302] V [ 132 , 386] ’  ) 

WRITE (8, 880) NW 

860  FORMAT  (  *  W(N’I1 * ) P [838 , 172] C [832 , 186] ’  ) 

WRITE (8 ,886) 

886  FORMAT  (  ’  P [846 , 172] Vll 1177771 1 1 177771 1 1 lOOO ’  ) 

WRITE (8, 870) NW 

870  FORMAT  (  ’  W (N ’ II ' ) P  [860 , 372] C [843 , 372]  ’  ) 

WRITE (8, 880) 

880  FORMAT  (  ’  P [867 , 372] VI 1 1 177771 1 1 lOO ’  ) 

ELSE 

WRITE (8, 800) NW 
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900  FORMAT  (  1  W  (N  ’  II  ’)  P  [  180 , 282]  V  [2-40 , 182]  ’  ) 

WRITE (0,905) 

906  FORMAT  (  ’  P [ 180 , 382] V [340 , 382] P [ 1 80 , 382] V [240 , 382] ’  ) 
WRITE (9, 910) NW 

910  FORMAT  (  ’  W (N ‘ II ’ ) P [1 13 , 272]  V [ 143 , 272]  '  ) 

WRITE (9,916) 

016  FORMAT  (  *  P  [138 , 268] V [] V [143 , 272] V[138,276]  ’  ) 

WRITE (0 , 920) NW 

020  FORMAT  (  *  W (N ’ II ’ ) P [690 , 172] V [816 , 172] ’  ) 

WRITE (9, 926) 

926  FORMAT  (  *  P [808 , 188] V [] V [816 , 172] V  [808 , 178]  ’  ) 

WRITE (0, 030) NW 

930  FORMAT  (  ’  W (N ’ I 1 ’ ) P [810 , 222] V [810 , 262] 1  ) 

WRITE (9, 936) 

036  FORMAT  (  ’  P [808 , 246] V [] V [810 , 252] V [614 , 246]  ’  ) 

WRITE (9 , 940) NW 

940  FORMAT  (  '  W(N ' 11 ’ ) P [810 , 372] V [580 , 372] ’  ) 

WRITE (9, 946) 

946  FORMAT  (  1  P [687, 388] V[]V [580, 372] V [687, 376] ’  ) 

WRITE (9 , 960) NW 

960  FORMAT  (  ’  W (N ’ II ’ ) P [160 , 372] V [128 , 372] V [128 , 390]  *  ) 
WRITE (9 ,966) 

966  FORMAT  (  ’  P  [124 , 381] V  [] V [128 , 300] V  [132 , 381 ]  ‘  ) 

WRITE (9 , 980) NW 

980  FORMAT  (  ’  W(N ’ II * )P [860, 172] C [843 , 172] ’  ) 

WRITE (9, 966) 

986  FORMAT  (  ’  P [867 , 172] Vl 11177771 1 llOO '  ) 

WRITE (9, 970) NW 

970  FORMAT  (  ’  W (N ’ II ’ )P [838 , 372] C [633 , 388] ’  ) 

WRITE (0,980) 

080  FORMAT  (  *  P [846, 372] Vl 11 177771 11 177770000’  ) 

END  IF 
RETURN 
END 

SUBROUTINE  REGISOUT 
C  SETS  VT126  BACK  INTO  TEXT  MODE 
LOGICAL  ALPHA 
DIMENSION  ALPHA (2) 

DATA  ALPHA  /  27,  »\*  / 

TYPE  lO,  ALPHA 
lO  FORMAT  (  ’+’,  6A1,  t  ) 

RETURN 

END 

SUBROUTINE  DELAY  (  NTICKS  ) 

TYPE  lOO, NTICKS 

lOO  FORMAT  (  ’ +S (T< ’ 13 ’ >) ’ ,  S  ) 

RETURN 

END 

SUBROUTINE  SCALEM  (  Y , IYABS , W , RANGE  ) 

DY  =  ((  W  -  Y  )  *  300.0  )  /  RANGE 
IDY  =  IFIX  (  DY  ) 

IYABS  =  IDY  +  90 

RETURN 

END 

SUBROUTINE  STIME  (  TIME , IXABS , RANGE , AXLEN , SHIFT  ) 
INTEGER  SHIFT 

DX  =  (  TIME  *  AXLEN  )  /  RANGE 
IDX  =  IFIX  (  DX  ) 

IXABS  =  IDX  +  SHIFT 

RETURN 

END 

SUBROUTINE  FRAME 

COMMON/SYST/DT , PSUP , POUT , TF , WBRL , VS , VO , IIN , TCYC 
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COMMON/NPRINT/NOUTl 
DIMENSION  T (lO) 

DELT  =  TF/IO.O 
DO  1=1,10 
T(I)=I*DELT 
END  DO 
TYPE  lO 

FORMAT  (  * !P[13O,30O]V[+8O,]V[,+6]  ’  ,  S  ) 

DO  1=1,0 
TYPE  12 

FORMAT  (  *  !P[,-6]V080,]V[,+6]  ’  ,  8  ) 

END  DO 
TYPE  20 

FORMAT  (  ’ !P[120,308]T"0.0" ’  ,  8  ) 

TYPE  21  ,  T (1) ,T(2) 

FORMAT  (  ’ !P[160,308]T" ’F6. 1 ' «P [220 , 308] T" ’F6. 1 ’ *  *  ,  8  ) 

TYPE  22  ,  T (3) , T (4) 

FORMAT  (  ’ !P[280,308]T" ’F5. 1 ’ "P [340 , 308] T* ’F6. 1 » *  *  ,  8  ) 

TYPE  23  ,  T (6) , T (8) 

FORMAT  (  ’ !P[400,308]T" ’F6. 1 ’ "P [480 , 308]  T" ’F6. 1 ’ " ’  ,  8  ) 

TYPE  24  ,  T (7) ,T(8) 

FORMAT  (  ’ !P[62O,308]T" ’FS. 1 ' "P [680 , 308] T" *F5. 1 ’ " '  ,  S  ) 

TYPE  26  ,  T (0) , T (lO) 

FORMAT  (  *  !P[84O,308]T" ’F5. 1 ’ "P [700, 308] T* ’F6. 1  * " *  ,  8  ) 
TYPE  30 

FORMAT  (  ' !P[360,430]T"TIME  (SEC.)*’.  *  ) 

TYPE  40 

FORMAT  (  '  !P[130,300]V[,-30] V[-6,]  '  ,  8  ) 

DO  1=1,0 
TYPE  60 

FORMAT  (  ’ !P[+6,]V[,-30]V[-6,] *  ,  8  ) 

END  DO 
TYPE  61 

FORMAT  (  ’ !P [70 , 180] T [+0 , +18] "OUTPUT" ’  ,  8  ) 

TYPE  63  ,  PSUP 

FORMAT  (  ’ !P [630 , 34] T [+0 , +0] "PSUP  =  "T* ’ F6 . 2 ’ »T"  PSIA"’,8) 
TYPE  64  ,  POUT 

FORMAT  (  1 !P[63O,62]T"P0UT  =  "T* ’F6 . 2 ’ "T"  PSIA"’,8) 

TYPE  66  ,  WBRL 

FORMAT  (  ’ !P[630,70]T"WBRL  =  ■T" ’ F6 . 2 ’ "T"  STD  LIT/MIN" ’ , $) 
IF (NOUTl . EQ . 1 ) THEN 
TYPE  80 

FORMAT  (  • !P  [200,6] T"*  OXYGEN  MOLE  FRACTION  ,  8  ) 

TYPE  82 

FORMAT  (  ’  !P[06,383]T"O.O,,P[O6,363]T''O.1"P[06,323]T''O.2"  ’8) 
TYPE  84 

FORMAT  (  ’  !P[06,203]T"0.3"P[86,283]T"0.4"P[e6,233]T"0.6"  ’  ,1 
TYPE  88 

FORMAT  (  ’ !P[06,203]T"0.8"P[06,173]T"0.7"P[06,143]T"0.8* ’ ,1 
TYPE  88 

FORMAT  (  ’ !P[06,113]T"0.0"P[06,83]T"1.0" ’ ,$) 

ELSE  IF  (NOUTl .  E(j .  2)  THEN 
TYPE  80 

FORMAT  (  ’ !P[20O,6]T"*  INLET  MASS  FLOWRATE  (KG/SEC)  **’,S) 
TYPE  82 

FORMAT  (  ’ !P[86,383]T"-.02"F [86 , 363] T" - . 01 *P [06 , 323] T"0 . O" ’ 
TYPE  84 

FORMAT  (  * !P[05,203]T" . Ol "P [06 , 283] T" . 02"P [06 , 233] T" .03" * ,8 
TfPE  88 

FORMAT  <  ’ !P[06,2O3]T" .O4"P[06, 173]T" .06"P[06, 143]T" .08" * ,8 
TYPE  88 

FORMAT  (  ’ !P[06,113]T".07"P[06,83]T".08"’,8) 

ELSE  IF (NOUTl .EQ. 3)  THEN 


TYPE  90 

90  FORMAT  (  ■ !P[290,6]T"*  OUTLET  MASS  FLOWRATE  (KG/SEC) 
TYPE  82 
TYPE  84 
TYPE  86 
TYPE  88 
ELSE 

GO  TO  67 
END  IF 
RETURN 
END 
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APPENDIX  C. 


These  are  the  corrected  equations  for  the  van  der  Viist  paper  (8).  These  are  for  the 
binary  mixture  at  total  pressure  (constant)  P.  The  six  unknowns  (A's  and  B's)  are  found  in 
terms  of  known  quantities  (G's  and  W’s).  See  the  paper  for  the  complete  details. 
Essentially,  six  equations  in  six  unknowns  are  solved.  The  equations  are: 


G0=  P  *  A, 

(1) 

Gj=  2  *  P  *A2  -A1+B1  +  2*P*B2  +  3*P2*B3 

(2) 

G.,=  3  *  P  *  A3  -  2  *  A,,  -  2  *  B2  -  3  *  2  *  P  *  B3 

(3) 

G3=  -  3  *  A3  +  3  *  B3 

(4) 

Wj  =  P  *  Aj  +  P2  *  A2  +  P3  *  A3 

(5) 

W2=  P  *  Bt  +  P2  *  B2  +  P3  *  B3 

(6) 

These  equations  give  the  following  solution: 


-6*W2.3*W1+2*G3*P3  +  3*GJ*P2  +  6*G,*P  +  9*a0 


■6*W2-6*W1+2*G!*P3  +  3*G2*P2  +  6*G|*Pt  12  *  G„ 
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